Estimating the seismic hazard in earthquake-prone regions, in order to assess the risk associated to nuclear facilities, must take into account a large number of uncertainties, and in particular our limited knowledge of the geology. And yet, we know that certain geological features can create site effects that considerably amplify earthquake ground motion. By combining the precision of physics-based simulations with the representativeness of deep neural networks, it becomes possible to quantify the influence of geological heterogeneity on the rendered earthquake space-time site response. This work demonstrates the advantages of adopting Factorized Fourier Neural Operator (F-FNO) to learn the Green’s operator of the elastodynamics partial differential equation, in order to quantify the influence of geological heterogeneity on the solution (the displacement or velocity field) at the free surface [1]. The trained F-FNO learns the relationship between 3D heterogeneous geologies and ground motions generated by the propagation of seismic waves through these geologies. The F-FNO is trained on the HEMEW-3D database, comprising 30000 high-fidelity numerical simulations of earthquake ground motion through generic geologies, performed by employing the high-performance code SEM3D [2, 3]. Next, a smaller database was built specifically for the Teil region (Ardèche, France), where a MW 4.9 moderate shallow earthquake occurred in November 2019 [4]. The F-FNO is then specialized on this database database with just 250 examples. Transfer learning improved the prediction error by 22 %. According to seismological Goodness-of-Fit (GoF) metrics, 91% of predictions have an excellent GoF for the phase (and 62% for the envelope). Ground motion intensity measurements are, on average, slightly underestimated.