Seismic hazard analyses in the area of a nuclear installation must account for a large number of uncertainties, including limited geological knowledge. It is known that some geological features can create site-effects that considerably amplify ground motion. Combining the accuracy of physics-based simulations with the expressivity of deep neural networks can help quantifying the influence of geological heterogeneities on surface ground motion.This work demonstrates the use of a Factorized Fourier Neural Operator (F-FNO) that learns the relationship between 3D heterogeneous geologies and time-dependent surface wavefields. The F-FNO was pretrained on the generic HEMEW-3D database with 30 000 samples. Then, a smaller database was built specifically for the region of the Le Teil earthquake (South-Eastern France) and the F-FNO was further trained with only 250 specific samples. Transfer learning improved the prediction error by 22 %.As quantified by the Goodness-Of-Fit (GOF) criteria, 90% of predictions had excellent phase GOF (62% for the enveloppe GOF). Although the intensity measures of surface ground motion were, in average, slightly underestimated by the FNO, considering a set of heterogeneous geologies always led to ground motion intensities larger than those obtained from a single homogeneous geology. These results suggest that neural operators are an efficient tool to quantify the range of ground motions a nuclear installation could face in the presence of geological uncertainties.