This dissertation is focused on creating the key ingredients of a numerical platform for reduced order modeling of soil structure interaction (SSI) problems. Specifically, a computer code is developed for forward simulation of wave propagation in the two-dimensional (plane-strain and axisymmetric) semi-infinite heterogeneous solid media. Perfectly matched layers (PMLs) are used for absorbing the outgoing waves. The computationally efficient symmetric hybrid PML formulation available for the plane-strain setting is extended to axisymmetric problems. The domain reduction method (DRM) is used for translation of the remote excitation within a PML-truncated medium. The methodologies are devised for using this finite element (FE) solver to (i) compute the soil impedance functions and (ii) the modified input motions (a.k.a. foundation input motion) of rigid and flexible interfaces embedded in heterogeneous half-spaces numerically. Existing semi-analytical solutions are used to verify these methods comprehensively. The verified framework is validated using data from a large-scale field test as well as centrifuge experiments. In order to demonstrate the framework's application: (i) the impedance functions and kinematic interaction transfer functions (KITFs) of a number of SSI problems---for which existing analytical solutions are limited---are computed; and (ii) the reduced order model of a buried structure in an elastic half-space is constructed. In order to avoid integro-differential equations, stable discrete-time filters are used for devising time-domain representations of the computed soil impedance matrix. The dynamic response of the resulting spatio-temporal reduced order model is compared against those obtained from solving the same problem using the direct modeling approach, and excellent agreement is observed.