In this work, we present a 2D elastodynamic model for the seismic response of subway tunnels embedded in a laterally inhomogeneous, multilayered geological region overlying the half-plane. To this end, a finite difference-boundary element methodology (FDM-BEM) is developed, with the latter method embedded in the former so as to capture near-site field effects. More specifically, the FDM is used for simulating in-plane elastic wave propagation from the underlying bedrock through the overlying soil deposits to the surface. A ‘box’ area is then defined within the original FDM mesh and contains lined tunnels. The ‘box’ is modeled by the BEM and its upper boundary coincides with the free surface of the geological deposit. This way, seismically-induced motions are imparted from the FDM mesh to the ‘box’ perimeter, so that the BEM may now be used to efficiently model the near-site layers which contain the tunnels. Verification studies are then successfully conducted for upward moving Gabor pulses, using the FDM alone, the present hybrid FDM-BEM and a hybrid FDM-finite element method formulation. Given that the FDM is defined in the time domain and the BEM in the frequency domain, the fast Fourier transform is used for linking these two constituent parts of the hybrid approach. This methodology is finally applied to a north–south geological cross-section of Thessaloniki, Greece, which contains two Metro tunnels placed directly below an important Roman-era monument known as the Arch of Galerius. Results are then given in the form of free-surface motions stemming from the Thessaloniki 5 July 1978 aftershock recorded at bedrock so as to establish the influence of the ongoing Metro line construction, now temporarily halted because of the economic crisis, on the free surface motions in the city centre where a number historical monuments besides the Arch of Galerius still survive.