In this work, we address the characterization of porosity and water saturation in a synthetic model of a shallow alluvial subsurface using frequency electromagnetic and seismic data. The inversion method employs a Gauss–Newton scheme, where the Jacobian of the merged seismic and electromagnetic data is formulated with respect to the spatially heterogeneous petrophysical parameters. This is made possible by introducing realistic petrophysical relationships, which significantly reduce the number of unknowns in the inverse problem and incorporate a strong prior correlation between the information contained in both data types regarding the subsurface composition. The results obtained show that this simultaneous joint petrophysical inversion produces reconstructions with clear improvements compared to independent petrophysical inversion. Indeed, it greatly enhances the spatial resolution of subsurface mapping, as well as the quantitative estimation of porosity and saturation.