\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRGEOS,Unicode,biblatex,published]{cedram}
\addbibresource{crgeos20250034.bib}
\newenvironment{Table}{\begin{table}}{\end{table}}
\newenvironment{noXML}{}{}
\def\defcitealias#1#2{}
\let\citetalias\textcite
\let\citepalias\parencite
\let\citep\parencite
\let\citet\textcite
\def\xcitealp#1#2{\citeauthor{#1}, \citelink{#1}{#2}}
\newcommand*{\citelink}[2]{\hyperlink{cite.\therefsection @#1}{#2}}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.95\linewidth}{#1}}
\def\botline{\\\hline}
\def\dollar{\$}
\def\tbody{\noalign{\relax}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\og{\guillemotleft}
\def\fg{\guillemotright}
\let\permil\textperthousand
\RequirePackage{etoolbox}
\usepackage{upgreek}
\usepackage{wasysym}
\usepackage{pifont}
\def\jobid{crgeos20250034}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\newcounter{runlevel}
\usepackage{multirow} 
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\def\back#1{}
\allowdisplaybreaks
\csdef{Seqnsplit}{\\}
\newcommand*{\hyperlinkcite}[1]{\hyper@link{cite}{cite.#1}}
\newenvironment{inftab}{}{}
\DOI{10.5802/crgeos.330}
\datereceived{2025-01-14}
\daterevised{2025-09-08}
\datererevised{2026-02-24}
\dateaccepted{2026-02-26}
\ItHasTeXPublished
}
{
\PassOptionsToPackage{authoryear}{natbib}
\documentclass[crgeos]{article}
\usepackage[T1]{fontenc} 
\usepackage{pifont}
\def\CDRdoi{10.5802/crgeos.330}
\def\href#1#2{\url[#1]{#2}}
 \def\citelink#1#2{\citeyear{#1}}
 \def\xcitealp#1#2{\citealp{#1}}
\makeatletter
\def\CDRsupplementaryTwotypes#1#2{}
 \def\selectlanguage#1{} 
 \def\hyperlinkcite#1#2{\citeyear{#1}}
\let\newline\break
}
\makeatother

\dateposted{2026-04-08}
\begin{document}

\begin{noXML}

\editornote{Article submitted by invitation}
\alteditornote{Article soumis sur invitation}

\CDRsetmeta{articletype}{review}

\TopicFR{Hydrologie, hydrog\'eologie}
\TopicEN{Hydrology, hydrogeology}
\TopicFR{G\'eophysique de surface}
\TopicEN{Surface geophysics}

\title{Review of integrating thermal, electrical, and seismic
geophysical methods for shallow groundwater modeling in the critical
zone}

\alttitle{Int\'egration des m\'ethodes g\'eophysiques thermiques,
\'electriques et sismiques pour la mod\'elisation des eaux souterraines
dans la zone critique}

\author{\firstname{Agn\`es} \lastname{Rivi\`ere}\CDRorcid{0000-0002-6002-3189}\IsCorresp}
\address{Mines Paris, PSL University, Centre for Geosciences and
Geoengineering, Fontainebleau, France}
\email[A. Rivi\`ere]{agnes.riviere@minesparis.psl.eu} 
\urladdr{https://agnes-riviere.github.io/}

\author{\firstname{Maxime} \lastname{Gautier}\CDRorcid{0000-0002-6416-1781}}
\addressSameAs{1}{Mines Paris, PSL University, Centre for Geosciences
and Geoengineering, Fontainebleau, France}
\email[M. Gautier]{maxime.gautier@minesparis.psl.eu}

\author{\firstname{Nicolas} \lastname{Radic}\CDRorcid{0009-0009-9275-294X}}
\addressSameAs{1}{Mines Paris, PSL University, Centre for Geosciences
and Geoengineering, Fontainebleau, France}
\email[N. Radic]{nicolas.radic@minesparis.psl.eu}

\author{\firstname{Ludovic} \lastname{Bodet}\CDRorcid{0000-0003-2271-3223}}
\address{Sorbonne Universit\'e, CNRS, EPHE, UMR 7619 METIS, 4 place
Jussieu, 75252 Paris 05, France}
\email[L. Bodet]{ludovic.bodet@sorbonne-universite.fr}

\author{\firstname{Alexandrine}\nobreakauthor\lastname{Gesret}\CDRorcid{0000-0002-6828-392X}}
\addressSameAs{1}{Mines Paris, PSL University, Centre for Geosciences
and Geoengineering, Fontainebleau, France}
\email[A. Gesret]{alexandrine.gesret@minesparis.psl.eu}

\author{\firstname{Roland} \lastname{Martin}\CDRorcid{0000-0003-3669-094X}}
\address{G\'eosciences Environnement Toulouse, Observatoire
Midi-Pyr\'en\'ees, Universit\'e de Toulouse, 14 Av. E. Belin, 31400,
Toulouse, France}
\email[R. Martin]{roland.martin@get.omp.eu}

\author{\firstname{Sylvain} \lastname{Pasquet}\CDRorcid{0000-0002-3625-9212}}
\addressSameAs{2}{Sorbonne Universit\'e, CNRS, EPHE, UMR 7619 METIS, 4
place Jussieu, 75252 Paris 05, France}
\email[S. Pasquet]{sylvain.pasquet@sorbonne-universite.fr}

\author{\firstname{Mathias} \lastname{Maillot}\CDRorcid{0009-0001-1908-510X}}
\address{Geological Survey of France -- BRGM, 3 Av. Claude Guillemin,
45100 Orl\'eans, France}
\email[M. Maillot]{m.maillot@brgm.fr}

\author{\firstname{Didier}\nobreakauthor\lastname{Renard}\CDRorcid{0000-0002-9402-0007}}
\addressSameAs{1}{Mines Paris, PSL University, Centre for Geosciences
and Geoengineering, Fontainebleau, France}
\email[D. Renard]{didier.renard@minesparis.psl.eu}

\author{\firstname{Pierre} \lastname{Weill}\CDRorcid{0000-0002-6201-1472}}
\address{Normandie Universit\'e, UNICAEN, UNIROUEN, CNRS, M2C, 14000
Caen, France}
\email[P. Weill]{pierre.weill@unicaen.fr}

\author{\firstname{Dominique} \lastname{Bruel}\CDRorcid{0000-0001-7873-8230}}
\addressSameAs{1}{Mines Paris, PSL University, Centre for Geosciences
and Geoengineering, Fontainebleau, France}
\email[D. Bruel]{dominique.bruel@minesparis.psl.eu}

\author{\firstname{R\'emi} \lastname{Valois}\CDRorcid{0000-0002-6447-7057}}
\address{Avignon Universit\'e, UMR EMMAH, Avignon, France}
\email[R. Valois]{remi.valois@univ-avignon.fr}

\author{\firstname{Jos\'e}\nobreakauthor\lastname{Cunha\nobreakauthor Teixeira}\CDRorcid{0000-0003-0305-4146}}
\addressSameAs{2}{Sorbonne Universit\'e, CNRS, EPHE, UMR 7619 METIS, 4
place Jussieu, 75252 Paris 05, France}
\email[J. Cunha Teixeira]{jose.teixeira@sorbonne-universite.fr}

\author{\firstname{Patrick} \lastname{Goblet}\CDRorcid{0000-0003-4540-4132}}
\addressSameAs{1}{Mines Paris, PSL University, Centre for Geosciences
and Geoengineering, Fontainebleau, France}
\email[P. Goblet]{patrickgoblet@free.fr}

\author{\firstname{Mark} \lastname{Noble}\CDRorcid{0000-0001-5638-6115}}
\addressSameAs{1}{Mines Paris, PSL University, Centre for Geosciences
and Geoengineering, Fontainebleau, France}
\email[M. Noble]{mark.noble@minesparis.psl.eu}

\keywords{\kwd{Hydrogeophysics}\kwd{Critical zone}\kwd{Transient
numerical models}\kwd{Inverse problem}\kwd{Uncertanty 
quantification}\kwd{Geostastistics}\kwd{Groundwater}}

\altkeywords{\kwd{Hydrog\'eophysique}\kwd{Zone critique}\kwd{Mod\`eles
num\'eriques transitoires}\kwd{Probl\`emes inverses}\kwd{Quantification
des incertitudes}\kwd{G\'eostatistique}\kwd{Aquif\`ere}}

\shortrunauthors

\begin{abstract} 
Groundwater systems are components of the Critical Zone, in dynamic
balance with the climate and human pressure. Water and heat fluxes
control biogeochemical processes and regulate both resource potential
and system vulnerability. Evaluating system responses to climate and
land use changes requires an estimation of subsurface structure, flow
dynamics, and the spatial distribution of recharge and discharge zones.

Classical hydrological data, such as water levels and river discharge,
are commonly coupled with unconventional methodologies, including heat
tracing, electrical resistivity tomography, and surface wave analysis.
Their coupling via time-lapse monitoring and inversion frameworks
enables the characterization of thermal, electrical, and mechanical
interactions, thus encouraging a transition from static structure
description to transient representations of groundwater.

We identify key challenges, including data limitations, model
uncertainty propagation, and the integration of transient variables and
parameters into inversion workflows. Petrophysical and geostatistical
approaches help overcome these issues by coupling geophysical and
groundwater models, quantifying spatio-temporal uncertainties, and
addressing scale change. Finally, choices regarding experimental
design, parameter reduction, dimensionality, model hypotheses, and
inversion scale must be assessed to balance parsimony with model
accuracy. These developments underscore the central role of
hydrogeophysics in advancing Critical Zone science and sustainable
groundwater management. Emphasizing transient processes is essential
for capturing how subsurface systems evolve over time in response to
environmental changes.
\vspace*{-3pt}
\end{abstract}

\begin{altabstract} 
Les aquif\`eres constituent une composante majeure de la Zone Critique,
en \'equilibre dynamique avec le climat et les pressions anthropiques.
Les flux d'eau et de chaleur contr\^{o}lent les processus
biog\'eochimiques et d\'eterminent la vuln\'erabilit\'e et le potentiel
de la ressource.

Aux donn\'ees hydrologiques classiques s'ajoutent des m\'ethodes
d'imagerie g\'eophysique permettant de visualiser l'eau dans le
sous-sol. Le tra\c{c}age thermique, la tomographie de r\'esistivit\'e
\'electrique et l'analyse d'ondes de surface fournissent des images
tridimensionnelles dont la r\'ep\'etition temporelle conduit \`a une
imagerie 4D, r\'ev\'elant la dynamique transitoire des aquif\`eres et
les variations de stockage d'eau et de chaleur.

Plusieurs d\'efis persistent : la raret\'e des donn\'ees, la
propagation des incertitudes et l'int\'egration des param\`etres
transitoires dans les sch\'emas d'inversion. Les approches
p\'etrophysiques et g\'eostatistiques contribuent \`a y r\'epondre en
couplant mod\`eles g\'eophysiques et hydrog\'eologiques, en quantifiant
les incertitudes et en traitant les changements d'\'echelle. Les choix
li\'es \`a la conception exp\'erimentale, \`a la r\'eduction des
param\`etres et aux hypoth\`eses de mod\'elisation doivent \'egalement
\^{e}tre \'evalu\'es pour assurer un compromis entre parcimonie et
fiabilit\'e.

Ces d\'eveloppements confirment le r\^{o}le central de
l'hydrog\'eophysique pour analyser l'\'evolution des aquif\`eres,
renforcer la surveillance des nappes et soutenir une gestion durable
compatible avec la transition \'energ\'etique. Ils contribuent
\'egalement \`a am\'eliorer la compr\'ehension des cycles
biog\'eochimiques et de l'\'evolution des structures et propri\'et\'es
de la Zone Critique.
\end{altabstract} 

\thanks{French Agence Nationale de la Recherche (ANR) under grant
ANR-23-CE01-0004 (project GWSBOUND), Institutes UMR 7619 METIS, UAR3455
ECCE TERRA, Sorbonne Universit\'e, CNRS, Mines Paris - PSL (PSL
University)}

%\input{CR-pagedemetas}

\maketitle

\vspace*{-2pt}

\twocolumngrid

\end{noXML}

\defcitealias{emerychangesupport2009}{ibid.}

\section{Introduction}\label{sec1}

\vspace*{-3pt}

The Earth's Critical Zone (CZ) is the near-surface layer where water,
soil, air, and living organisms interact to support terrestrial
ecosystems and human society. Among its components, groundwater (GW)
represents the largest reservoir of freshwater
\citep{singha_importance_2022,lee_meanings_2023}. It plays a
fundamental role in Earth's hydrogeochemical cycles by redistributing
water, energy, solutes, and contaminants throughout the subsurface
\citep{brunke_ecological_1997,stegen_influences_2018,gomez-gener_towards_2021}. 
In addition to sustaining ecosystems and surface water (SW) bodies
critical to biodiversity, GW is a vital resource for global water,
energy, and food security
\citep{de_amorim_nexus_2018,castilla-rho_social_2017,sacco_groundwater_2024}.

Despite its significance within the Critical Zone (CZ), GW is
frequently overlooked due to the difficulty and high cost of subsurface
access, and its characteristically slow, diffuse dynamics
\citep{chavez_garcia_silva_multi-decadal_2024}. {Hydrogeophysical
methods} provide efficient, non-intrusive tools to estimate GW physical
properties in time and space. However, key challenges in CZ
hydrogeophysics include: (i) characterizing the structure and behavior
of GW systems and their interactions with other CZ compartments, (ii)
quantifying the temporal and spatial distribution of water, energy, and
nutrients or contaminants, and (iii)~predicting GW responses to
climatic variability and anthropogenic influences.

Shallow GW systems are influenced by both climatic factors
\citep{cuthbert_global_2019, fan_global_2013,
condon_evapotranspiration_2020} and anthropogenic pressures such as
water abstraction \citep{dixon_subsidence_2006,
rahardjo_effects_2010,waltham_sinkholes_2004}, geothermal exploitation
\citep{bidarmaghz_impacts_2021,menberg_opportunities_2025,bayer_geothermal_2019}, 
and contaminant discharge \citep{stigter_groundwater_2023}. In this
context, understanding long-term trends and the spatial dynamics of GW
recharge, SW--GW exchanges, flow paths, residence times, and associated
thermal regimes is essential for anticipating future conditions and
managing current impacts \citep{maxwell_interdependence_2008,
de_marsily_quantitative_1986, enemark_hydrogeological_2019,
boulton_rivers_2006,
singha_importance_2022,gleeson_global_2020,fan_global_2013}. 
Temperature is increasingly recognized as a key tracer in shallow GW
systems, providing critical insights into water quality, recharge
dynamics, and geothermal potential
\citep{neidhardt_impact_2023,xu_effect_2021,benz_global_2024}. Although
only a limited number of studies have focused on temperature evolution
within GW, recent investigations underscore its sensitivity to climatic
variability and anthropogenic pressures \citep{benz_global_2024}. Water
and temperature are fundamental vectors of transport and transformation
within the CZ, shaping biogeochemical processes that influence CZ
functioning \citep{safieddine_surface_2025}, altering water quality
through complex biogeochemical feedbacks \citep{neidhardt_impact_2023},
and affecting stream and GW biodiversity \citep{land_groundwater_2023}.

Unraveling the four-dimensional (4D) distribution of water and heat
fluxes at the GW boundaries (i.e.\ 3D spatial ${+}$\ temporal), particularly
at the Vadose Zone (VZ) and SW--GW interfaces
\citep{hermans_advancing_2023} remains a major challenge for accurately
assessing the components of the water and energy balances and achieving
water budget closure, particularly concerning climate changes and land
use \citep{bloschl_twenty-three_2019,condon_global_2021,
xie_uncertainty_2018,xie_uncertainty_2019,levin_uncertainties_2023}. VZ
regulates recharge, vertical water movement, and the exchange of energy
and solutes, while SW--GW interactions control hydrological connectivity
as well as the magnitude and directionality of both water and thermal
exchanges. These efforts are further complicated by the non-linear,
transient nature of GW recharge processes, including capillary rise
across the water table \citep{vereecken_value_2008}, shifting GW
divides that can induce lateral flows between riverbanks
\citep{texier_numerical_2022}, and temporally variable SW--GW
connectivity. SW--GW can fluctuate between gaining, losing,
transitional, and disconnected states, introducing uncertainty into
water and energy budget estimation \citep{brunner_advances_2017,
brunner_disconnected_2011, riviere_experimental_2014}
(Figure~\ref{fig:Gw_bound}).

\begin{figure*}
\includegraphics{fig01}
\vspace*{2pt}
\caption{\label{fig:Gw_bound}(a) Watershed-scale view illustrating GW
flow paths, geological structure, SW--GW interactions, and human
influences, with idealized GW divides and flowpaths 
\citep[after][]{toth_theoretical_1963}. (b) Hillslope-scale 3D
view with subsurface hydrofacies heterogeneity meshes, model parameters
(intrinsic permeability $k$, porosity $n$, specific yield $Sy$, water
retention curve WRC
\citep{van_genuchten_closed-form_1980,mualem_new_1976,brooks_properties_1966,fredlund_unsaturated_2006}, 
thermal conductivity $\lambda$, the density of the porous media $\rho$
and heat capacity $C$). Parameters may vary continuously or discretely
across the domain, reflecting natural heterogeneity and modeling
assumptions. (c) Output of advanced GW simulation following data
integration and inversion (Inv), revealing realistic SW--GW flow paths,
residence time distribution, and estimates of water and heat budgets.
After inversion, spatial uncertainty in model parameters and flow
directions can be quantitatively assessed, reflecting the limitations
imposed by sparse or uneven field data. FWD: forward models; Inv:
inversion approaches.}
\vspace*{2pt}
\end{figure*}

GW flow systems remain poorly understood due to their hidden subsurface
nature, often described as a terra incognita
\citep{kleinhans_terra_2005, mcdonnell_beyond_2017,
grant_frontier_2017} (Figure~\ref{fig:Gw_bound}). Shallow GW systems
exhibit a rich diversity in the composition and topology of internal
structures, shaped by interactions among geological, hydrological, and
climatic processes (see Appendix~\ref{app:geology}). These systems are
underlain by regolith, which can be categorized into two types:
allochthonous regolith, composed of unconsolidated sediments
transported through geomorphological mechanisms, and autochthonous
regolith, formed through in-situ weathering and transformation of the
underlying bedrock. The evolution of regolith and bedrock varies
significantly across spatial and temporal scales, complicating the
accurate representation of flow and transport dynamics (further
detailed in Appendix~\ref{app:geology}). Despite the advancement of
hydrogeological modeling, GW management, including well regulation,
streamflow depletion, recharge enhancement, and development of
geothermal resources, often operates under significant data
constraints. These limitations require simplifying assumptions that can
obscure spatial heterogeneity, misrepresent boundary conditions, and
ultimately lead to non-unique or unreliable model \mbox{interpretations}
\citep{casillas-trasvina_bayesian_2024}. However, the hydrological
sciences have historically prioritized hydraulics and mathematical
modeling over geological realism, frequently yielding precise solutions
to ill-posed problems. Inadequate conceptual geological models and poor
estimations of aquifer parameters still contribute to major sources of
uncertainty in flow modeling \citep{rojas_assessment_2010,
refsgaard_review_2012, enemark_hydrogeological_2019,
yin_accounting_2021, tran_influence_2025}.\looseness=1

While many components of the hydrological cycle, such as river
discharge, precipitation, piezometric levels, and even GW pumping
(though often not routinely monitored), can be measured with reasonable
accuracy, fluxes at GW boundaries, including recharge
\citep{rawls_infiltration_1993,halloran_heat_2016,zheng_response_2019,vereecken_soil_2022} 
and SW--GW exchanges
\citep{fleckenstein_groundwater-surface_2010,
krause_inter-disciplinary_2011,flipo_continental_2014,
barthel_hess_2014,boano_hyporheic_2014,irvine_groundwater-surface_2024}, 
cannot be directly observed at the hillslope scale. Although methods
such as seepage meters \citep{rosenberry_history_2020}, lysimeters
\citep{saaltink_analysis_2020,balugani_lysimeter_2023,saaltink_analysis_2020}, 
single-well tracer tests \citep{paradis_singlewell_2022}, and active
Distributed Temperature Sensing (DTS) experiments
\citep{chang_application_2024,briggs_comparison_2012} are often
described as ``direct measurements'' of these fluxes
\citep{rosenberry2008field}, they actually infer exchange rates by
inverting data. Each inversion requires modeling assumptions and is
sensitive to site-specific conditions. While informative at the local
scale, such measurements are rarely representative of larger domains.

At the hillslope or watershed scale, SW--GW exchanges are often
depicted using a one-dimensional conductance term in watershed models.
This representation neglects heterogeneity and assumes steady-state
conditions and places the GW divide at the midpoint of the river
\citep{gianni_rapid_2016}. This conceptual limitation is illustrated in
Figure~\ref{fig:Gw_bound}a, which shows how GW flow paths and divides
are idealized without accounting for hydrofacies variability. Notably,
if the transient state is not taken into account, predictions of
exchange fluxes will be unreliable, especially during extreme events
\citep{tripathi_review_2021}. This limitation impacts water fluxes as
well as heat and solute transport processes
\citep{lemoubou_effects_2023}. In reality, GW head 
\mbox{distributions} follow
physical laws that yield complex spatial patterns, which are not easily
reproduced with analytical solutions or simplified geometries.
Consequently, assigning boundary conditions based on topography or
interpolation of piezometric level carries a significant risk of
physical inconsistency \citep{de_marsily_dealing_2005}. A major
unresolved issue remains the inherent incompleteness of the data
required to constrain these models accurately.\looseness=1

Figure~\ref{fig:Gw_bound} provides a conceptual overview of GW systems
at multiple scales, illustrating (a) watershed-scale flow paths and
human influences, (b) hillslope-scale subsurface heterogeneity and
hydrogeological parameters, and (c) the output of advanced GW
simulations and inversion approaches, which estimate GW flow paths,
residence times, SW--GW exchanges, and recharge across the domain. The
fundamental challenge lies in solving flow and transport equations
within a complex geological system that remains only partially
characterized \citep{anderson_applied_2015}. Parameters may vary
continuously or discretely across the domain. In both cases, the
subsurface domain is discretized into a mesh of cells in the model.
Each mesh assigned a set of parameters that defines the hydraulic and
thermal behavior of the medium (Figure~\ref{fig:Gw_bound}b). These
parameters such as the saturated hydraulic conductivity (or
permeability), porosity, the relative permeability function, water
retention curve, unsaturated permeability function, specific heat (heat
capacity/density), thermal conductivity, and density of every phase
(i.e., solid particles, water, and air) are typically linked to a
specific hydrofacies, which represent geologically coherent units with
distinct properties. While this framework enables spatial variability
in model inputs, it also depends heavily on how well the hydrofacies
are conceptualized, parameterized, and distributed
\citep{zinn_when_2003,zhan_subsurface_2023}, underscoring the
importance of robust geological characterization. Although
heterogeneity is incorporated into models at a certain scale, it is not
directly observed or mapped in detail at the scale of the cell of the
model (Figure~\ref{fig:Gw_bound}b). Instead, its spatial distribution
is typically inferred from sparse measurements through inverse modeling
or interpolated through a geostatistical framework, introducing further
uncertainty into the representation of subsurface complexity
\citep{condon_global_2021, zhou_inverse_2014}.

As a result, GW simulations often reflect a compromise between
\citep{gupta_towards_2012}: (i) available data and its spatiotemporal
resolution \citep{irvine_groundwater-surface_2024, condon_global_2021,
de_marsily_comment_1992, konikow_ground-water_1992,
andreassian_hydrogeological_2023, carrera_discussion_1993}; (ii) the
complexity of physical processes; (iii)~model capacity to represent
heterogeneity and boundary conditions; and (iv) the inclusion of
uncertainties in model inputs, structure, parameters, and observations
used to constrain model parameter sets \citep{condon_global_2021,
zhou_inverse_2014, sagar_direct_1975, tarantola_popper_2006}
(Figure~\ref{fig:Gw_bound}c). 

Most hydrological models are calibrated using two conventional
observation types: hydraulic heads and river discharge. However, it is
established that these conventional observations alone are insufficient
to accurately parameterize subsurface heterogeneity and predict flow
paths and transport dynamics, especially in geologically complex
settings and near the SW \citep{de_marsily_dealing_2005}. Calibration
from hydraulic heads alone often leads to non-unique solutions
\citep{beven_undermining_2006, doherty_modeling_2011}. This issue
arises because distinct hydrofacies can produce similar hydraulic
responses but different flow paths.

To illustrate the hydrogeological consequences of sedimentary
heterogeneity, Figure~\ref{fig:sed_setup} presents a hydrogeological
simulation of GW flow paths and hydraulic gradients within a synthetic
geological model of a meandering fluvial system. The geological model,
developed using the \href{https://flumy.minesparis.psl.eu/}{Flumy}
process-based sedimentary simulator, mimics the deposits in a
simplified alluvial plain
\citep{weill_process-based_2013,bhavsar_stable_2024}. The facies
simulated by Flumy were implemented in a GW flow model to assess the
impact of mud plugs, which are low-permeability deposits formed in
abandoned meanders, on subsurface flow dynamics. Modeling details are
provided in Appendix~\ref{app:geology}. The results show that while the
overall hydraulic head distributions in the homogeneous and
heterogeneous models appear nearly identical, the presence of mud plugs
significantly alters GW flow directions and locally doubles the
velocity beneath these features. These insights underscore the
limitations of relying solely on hydraulic head data for model
calibration.

\begin{figure*}
\includegraphics{fig02}
\caption{\label{fig:sed_setup}Sedimentary and hydrogeological setup
(facies and hydrological boundary conditions): (a) homogeneous
sedimentary model, (b) heterogeneous sedimentary model, (c,~d)
zoom of the simulated hydraulic head (black lines) and Darcy flow
velocity fields of both models, (e,~f) slice view below the
mudplug at $t=15$ days. The map-view locations of profiles A--B (panels
a and c) and A$'$--B$'$ (panels b and d) correspond to the vertical
slices shown in panels (e) and (f).}
\end{figure*}

Conventional observations coupled with unconventional data (heat or
biogeochemical tracer, geophysical data) allow addressing these
limitations of GW modeling \citep{de_marsily_dealing_2005,
schilling_beyond_2019,casillas-trasvina_bayesian_2024}. Concentration
and temperature data are more commonly available and jointly used with
the hydraulic head in inversion processes \citep{ellison_can_2025}.
Recent reviews demonstrate the role of geophysics in hydrological
studies, for characterizing underground hydrofacies, constraining the
GW architecture, and estimating key hydrological properties such as
hydrodynamic and thermal parameters
\citep{hubbard_hydrogeological_2000, parsekian_multiscale_2015,
binley_dc_2005, hermans_advancing_2023, wagner_overview_2021,
dumont_geophysics_2024}. These methods provide spatially extensive,
non-invasive insights into subsurface architecture and hydrofacies
distributions \citep{hubbard_hydrogeological_2000,
parsekian_multiscale_2015, binley_dc_2005, hermans_advancing_2023}.
\citet{beaujean_calibration_2014}, \citet{copty_testing_2016}, and
\citet{irving_stochastic_2010} show that integration of geophysics 
results inside GW models can significantly improve the accuracy of the
model compared to using hydrogeological data alone. Common geophysical
approaches such as electrical resistivity tomography (ERT), seismic
refraction tomography (SRT), and temperature tracer are widely applied
in transient states and complex, heterogeneous geological environments
at the hillslope scale.

Understanding and modeling shallow groundwater systems in heterogeneous
environments remains a major challenge in hydrology and CZ science.
This review explores how the integration of geophysical methods,
specifically electrical, thermal, and seismic techniques, can enhance
our ability to characterize and model these complex systems. We
synthesize recent advances, highlight methodological developments, and
propose directions for future research at the interface of geophysics
and hydrology.

We begin by outlining the fundamental principles and applications of
these three geophysical methodologies, highlighting their respective
strengths and limitations in the context of heterogeneous shallow GW
investigations. A particular focus is placed on the potential of these
methods to reveal 4D GW fluxes, which are often poorly constrained in
heterogeneous shallow GW. Recent advances in incorporating geophysical
data into groundwater models are then discussed, with a particular
focus on the development of various inversion approaches that optimize
data assimilation and improve model reliability.

Finally, we identify emerging methodological developments and propose
key research directions for integrating hydrogeophysical data into GW
models, thereby bridging the gap between theory and practice. These
include strategies for effective field data monitoring, the development
of scalable hydrogeophysical models, and robust approaches for
quantifying uncertainty. 

\section{Geophysical methods, petrophysics, and geostatistics}
\label{sec2}

To identify hydrofacies parameters, common hydrological practices
include laboratory measurements, pumping tests
\citep{theis_relation_1935, papadopulos_drawdown_1967,
jacob_nonsteady_1952,black_determination_1981}, GW time-series
analyses, such as data in the spectral domain
\citep{gelhar_groundwater_1974, guillaumot_frequency_2022,
pedretti_scale_2016, schuite_improving_2019}, or responses to Earth or
atmospheric tides \citep{klonne_periodischen_1880,
rau_technical_2020,meinzer_ground_1939, merritt_estimating_2004,
mcmillan_utilizing_2019, thomas_hydro-mechanical_2024,
valois_technical_2024}. Critics argue that these methods are inherently
local, dependent on the design of pumping test experiments, piezometer
placement, and site-specific conditions. There is ongoing debate about
the spatial representativeness and upscaling of these parameters from
local measurements to hydrogeological models at the hillslope scale.
However, these data, combined with detailed geological descriptions,
provide valuable constraints for hydrogeological models, improving the
calibration and validation of model parameters. These data are commonly
used as prior information in hydrogeological modeling, supporting the
range of the physical parameter values of the subsurface. Rocks and
unconsolidated materials are often described as a mix of various
components: rock matrix and multi-fluids (air, water). These
characteristics control the physical properties such as the hydraulics,
thermal, electrical, and seismic parameters. Therefore, geophysical
methods can assist in characterizing hydrofacies through direct
modelling and inverse problems solving \citep{rubin_stochastic_2005}.

\subsection{Thermal, electrical, and seismic methods}

Equations governing the forward model (FWD) of each method are given in
the Appendix~\ref{app:equations}.

\subsubsection{Heat as a flow tracer}

While temperature is widely acknowledged as a significant state
variable with extensive hydrological and hydrochemical impacts, the
past decade has seen an increased recognition of the value of
temperature measurements to trace subsurface flow paths, especially for
GW flow, vadose zone and SW--GW interaction
\citep{anderson_heat_2005,kurylyk_theory_2019,halloran_heat_2016}.
Owing to the strong coupling between flow and heat transport,
temperature has long been used as a natural tracer to infer subsurface
flow processes, beginning with foundational studies in the 1960s
\citep{suzuki_percolation_1960, stallman_steady_1965,
bredehoeft_rates_1965}. In such systems, the advective term in the heat
transport equation is directly proportional to specific discharge
(Appendices~\ref{app:flow_eq}, \ref{app:heat_eq}). The theoretical
basis stems from classical heat conduction work by
\citet{carlaw_conduction_1959} and has been extended to account for
advective heat transport under transient conditions
\citep[e.g.][]{campbell_probe_1991, heitman_field_2003,
briggs_comparison_2012, simon_numerical_2021}. The increasing
availability of affordable temperature sensors facilitates
high-resolution monitoring of thermal dynamics in rivers, soil, and
aquifers, which provides a natural signal for tracing water movement
(Figure~\ref{fig:thermic_method}a,~b). Temperature fluctuations,
whether diurnal, seasonal, or induced by artificial sources, can be
leveraged to infer SW--GW interactions, quantify the GW recharge,
delineate subsurface flow paths, and quantify geothermal potential
\citep{hatch_quantifying_2006, keery_temporal_2007,
kurylyk_theory_2019, saar_review_2011, anderson_heat_2005,
constantz_heat_2008, rau_analytical_2010, healy_using_2002}.

\begin{figure*}
\vspace*{4pt}
\includegraphics{fig03}
\vspace*{4pt}
\caption{\label{fig:thermic_method}Using heat as a tracer: (a)
Temperature time series recorded in the Orgeval CZ observatory (France)
\citep{tallec_observatoire_2015,riviere_revue_2018} for aquifer (black
line), bank (gray line), and stream (light-grey line) over a decade. 
(b) Annual temperature regime plot: light blue and red points indicate
the daily minimum and maximum temperatures, respectively; the
gray-shaded region represents the range between the 0.05 and 0.95
quantiles, and the dark blue line shows the mean. (c) Simulation of a
temperature profile as a function of depth for a homogeneous soil, and 
(d) a heterogeneous soil, for aquifer discharge (orange line) and
aquifer recharge (blue line). (c) and (d) showing deviations from the
geothermal gradient caused by surface zone and convection in the
geothermal zone. Aquifer recharge (downward movement) results in
concave upward profiles, whereas aquifer discharge (upward movement)
results in convex upward profiles.}
\vspace*{4pt}
\end{figure*}

The use of temperature as a hydrological tracer relies on two main
strategies: passive and active thermal tracing. Passive approaches take
advantage of natural temperature fluctuations, typically diurnal or
seasonal surface signals, and analyze their propagation through the
subsurface to infer GW flow characteristics
(Figure~\ref{fig:thermic_method}a,~b). Signal attenuation and phase
shifts with depth (Figure~\ref{fig:thermic_method}a) allow estimation of
vertical fluxes and thermal diffusivity, particularly in streambeds and
in the VZ \citep[e.g.][]{hatch_quantifying_2006,
keery_temporal_2007, kurylyk_theory_2019, constantz_heat_2008}. Heat is
a particularly effective tracer for revealing subsurface
heterogeneities, as its transport is highly sensitive to variations in
thermal properties and flow paths. 

Figure~\ref{fig:thermic_method}c and d present transient temperature
profiles produced with the Ginette model
\citep{riviere_experimental_2014,riviere_experimental_2019,
riviere_ginette_2020,grenier_groundwater_2018}, 
a coupled water- and heat-diffusivity equations, for both homogeneous
and two-layer domains, demonstrating the influence of \mbox{subsurface}
\mbox{heterogeneity} on heat transfer. {During aquifer recharge}, surface
water infiltrates downward, and advective heat transport produces a
concave-upward temperature profile characteristic of downward flow. The
infiltrating water may be warmer in summer or cooler in winter, but the
curvature remains diagnostic of recharge, while discharge creates a
convex or compressed profile near the surface
(Figure~\ref{fig:thermic_method}c and d) \citep{riviere_thermal_2020}.
These advective processes stand in clear contrast to purely conductive
heat transfer, where temperatures vary only by diffusion and follow a
nearly linear trend with depth. In the heterogeneous model, the
presence of a shallow, permeable layer accelerates advective transport
during recharge and discharge events. As a result, the thermal profile
diverges from the expected recharge pattern and instead bends
convex-upward. \citet{kurylyk_heat_2017} demonstrates that heat can
serve as an effective hydrologic tracer in both shallow and deep
heterogeneous media, offering analytical tools based on the work of
\citet{shan_analytical_2004} and field data to improve subsurface flow
characterization. Thermal tracing is widely used to assess SW--GW
interactions, as reviewed by \citet{anderson_heat_2005},
\citet{saar_review_2011}, and \citet{kurylyk_theory_2019}. However, its
accuracy relies on having enough advective heat transfer. Recommended
intrinsic permeability values range from $10^{-13}$ to $10^{-10}$
m\tsup{2} \citep{cucchi_estimating_2021, smith_thermal_1983,
riviere_thermal_2020}.\break If the conductivity falls below this range, it
is not possible to infer intrinsic permeability from thermal data, as
heat transport is dominated by conduction.\looseness=1

In the VZ, temperature data have also been used to estimate recharge
rates \citep[e.g.][]{tabbagh_determination_1999, benderitter_flow_1993,
bechkit_monitoring_2014, halloran_heat_2016}. However, inversion under
unsaturated conditions requires a fully coupled thermo-hydraulic
framework, involving soil-specific relationships for thermal
conductivity and water retention
\citep[e.g.][]{van_genuchten_closed-form_1980, mualem_new_1976,
brooks_properties_1966, fredlund_unsaturated_2006}. However, no single
model yet exists for thermal conductivity in VZ
\citep{cosenza_relationship_2003}. Thus, testing and comparison with
other methods under a variety of conditions, such as low saturation,
high saturation, and active infiltration at various rates, is an
obvious next step to improve the VZ heat tracer methodology
\citep{halloran_heat_2016}. Quantifying vertical GW discharge fluxes
necessitates measurement approaches with spatial resolutions consistent
with the scale of subsurface heterogeneity and the gradients driving
vertical flow \citep{tabbagh_determination_1999}.

Recent advances in measurement technologies have significantly expanded
the use of thermal methods to characterize GW flow and heterogeneity.
Fiber-optic distributed temperature sensing (FO-DTS) enables
high-frequency, high-resolution temperature measurements over long
spatial domains \citep[e.g.][]{tyler_environmental_2009,
hermans_geophysical_2014}. Applications include resolving small-scale
hydrofacies variability and preferential flow paths in alluvial
systems. In active FO-DTS, heat is applied along the fiber, and
temperature changes are used to infer thermal and hydraulic properties
through the inversion. This method has proven particularly effective in
high-flux environments, where passive signals may be obscured.
\citet{briggs_comparison_2012} and \citet{simon_numerical_2021}
demonstrated how thermal signals correlate with known hydrofacies and
enhance resolution of vertical water fluxes in the context of SW--GW
exchanges.

Drone-based thermal infrared imaging provides complementary surface
temperature data across broad spatial extents. In Arctic catchments and
braided river systems, UAV-mounted TIR cameras have been used to detect
GW discharge zones and identify ecologically significant thermal
habitats \citep[e.g.][]{hare_continental-scale_2021}. These systems
have been validated against FO-DTS and in-situ sensors, providing an
effective platform for regional thermal mapping.\looseness=1

Active thermal tracing offers a controlled approach to investigating
subsurface flow by introducing a heat source, via heated probes,
injected water, or fiber-optic cables {surrounded by a wire acting
as a resistance that can heat the medium}, and monitoring its spatial
and temporal propagation. This approach is advantageous where natural
thermal gradients are insufficient to resolve flow dynamics. Moreover,
imposing a well-defined heat source significantly reduces uncertainties
related to external forcing and simplifies the inversion process. For
instance, studies using actively heated DTS have demonstrated the
ability to estimate SW--GW exchanges \citep{briggs_actively_2016,
simon_combining_2022}. Analytical solutions, such as the moving
instantaneous line source model, can then be applied to derive thermal
conductivity and SW--GW fluxes with relatively low uncertainty
\citep{simon_numerical_2021, simon_combining_2022}. Similarly, thermal
response tests using actively heated fiber-optic cables in boreholes
have led to precise estimation of GW flow rates and the
characterization of borehole properties by modeling heat transfer in a
controlled environment \citep{zhang_estimation_2023}.
\citet{chang_application_2024} applied this approach for characterizing
subsurface heterogeneity and flow dynamics in a tidal-influenced
coastal aquifer. The thermal tracer needs to be used for a long period
of time, depending on the thermal delay effect
\citep{palmer_thermal_1992}. Passive temperature sensing might be
better for long-term and wide-ranging monitoring
\citep{furlanetto_fiber_2024}.

Thermal methods offer distinct advantages for studying flow in
heterogeneous media (Figure~\ref{fig:thermic_method}c,~d). Unlike
hydraulic head or solute concentration data, temperature signals can
directly characterize the direction, magnitude, and timing of flow
across hydrofacies boundaries. However, the nature of heterogeneity,
involving abrupt changes in lithology, porosity, hydraulic
conductivity, thermal conductivity, and water retention properties,
poses different challenges when interpreting heat transport in
subsurface systems. Hydraulic conductivity can vary dramatically, which
complicates the estimation of GW flow rates based solely on temperature
\mbox{measurements} \citep{constantz_heat_2008}. This necessitates the
integration of temperature observations with hydraulic gradient
measurements for robust flow characterization
\citep{cucchi_estimating_2021}.

Joint heat and solute tracer experiments have been proposed to reduce
uncertainties in transport predictions by exploiting their
complementary behaviors, offering a more comprehensive understanding of
subsurface heterogeneity and improving model accuracy
\citep{hoffmann_heterogeneity_2019}. However, complex three-dimensional
flow fields can still introduce significant errors in the estimation of
thermal fluxes and effective thermal diffusivity, underlining the need
for modeling approaches that explicitly consider spatial variability
and flow directionality \citep{sommer_impact_2013,reeves_impacts_2016}.

\subsubsection{Electrical resistivity tomography}

The electrical conductivity (or resistivity) of geological materials is
governed by the lithological properties of the rock matrix, as well as
the chemistry, volume, and temperature of internal fluids
\citep{hayley_low_2007}. Consequently, electrical investigation methods
provide a direct window into these subsurface parameters. ERT is a
preferred near-surface geophysical technique, valued for its
non-invasive nature and ease of implementation
\citep{everett_near-surface_2013}. Its widespread adoption is further
supported by various robust data processing software packages (see
Appendix~\ref{app:ert}) \citep{loke_rapid_1996, blanchy_resipy_2020,
rucker_pygimli_2017}. As illustrated in Figure~\ref{fig:ert_method},
the typical ERT workflow follows three stages: (i) installing electrodes
at the ground surface or within boreholes, (ii) collecting apparent
resistivity ($\rho_{a}$) data using a variety of quadrupole
configurations, and (iii) solving an inverse problem to generate a final
electrical resistivity model of the subsurface.

\begin{figure*}
\includegraphics{fig04}
\caption{\label{fig:ert_method}ERT method consisting of apparent
resistivity measurements during a field survey, data processing or
inversion, and interpretation phase. Streamlines and depth of
investigation (blue circle) for a given quadrupole \textit{ABMN} are
shown.}
\end{figure*}

ERT consists in inducing an electrical current ($I$) in the ground
using two injection electrodes (Figure~\ref{fig:ert_method}). The
resulting voltage difference ($\Delta{}U$) is measured by two potential
electrodes. These four electrodes form a quadrupole. The apparent
resistivity $\rho_{a}$ is then derived from the current, the voltage
difference, and the geometrical factor, which depends on the
arrangement of the quadrupole's electrodes. The apparent resistivity
represents an average value of the subsurface, assuming a homogeneous
half-space. {Increasing the spacing between the electrodes allows
sounding deeper, but the resulting value remains an integration over
the entire volume of ground between the injection electrodes.} All
quadrupoles of a profile are characterized by a theoretical depth of
investigation defined by \citet{roy_depth_1971} and
\citet{barker_depth_1989}. A 2D distribution of $\rho_{a}$ along a
transect can be obtained by combining several combinations of
``single'' measurements.

Measured apparent resistivity values are then processed and inverted to
eventually obtain a resistivity model. This model can be interpreted in
terms of subsurface structure \citep{guerin_geophysical_2009,
boura_improving_2024}, lithology
\citep{guerin_geophysical_2004,bentley_two_2004}, water content,
\citep{tabbagh_measurement_2002,tabbagh_determination_1999,
bai_groundwater_2021}, hydraulic conductivity \citep{rehfeldt_field_1992}, 
salt content \citep{hayley_time-lapse_2009}, transport parameters of
tracer plumes \citep{lekmine_quantification_2017,
benderitter_heat_1982, giordano_time-lapse_2017,
hermans_geophysical_2014}, or water level
\citep{bentley_accuracy_2002,wu_application_2023}. Inferring parameters
from geological and hydrofacies features is possible using
petrophysical models from the literature or site-specific information
(Section~\ref{sec:petrophy}). For dynamic parameters or variables, a
time-lapse approach is required (cf.\ Section~\ref{sec:transient}). Data
errors significantly impact both the inversion stability and the
subsequent interpretation of the resulting resistivity images
\citep{tso_improved_2017, tso_field_2019}. Often, field data errors are
estimated by stacked or reciprocal measurements performed with the
resistivity meter \citep{tso_improved_2017, tso_field_2019,
dahlin_short_2000, wilkinson_practical_2012}. Error estimation is
essential in the inversion process because underestimating error levels
can result in model overfitting, while overestimating them can lead to
underfitting, both of which complicate the interpretation of results.
ERT output models are particularly challenging to interpret due to the
unconstrained nature of ERT images, which are affected by survey
configuration, field conditions, subsurface heterogeneity, and the
inherent non-uniqueness of the inversion solution. Some tools, such as
the PyMERRY code \citep{gautier_pymerry_2024}, can be employed to
quantify uncertainties in ERT images in addition to classical
sensitivity computation \citep{loke_tutorial_2025}.\looseness=1

\subsubsection{Active seismic methods}

Seismic approaches are widely employed at various scales in
hydrogeology; however, they primarily focus on delineating subsurface
contrasts rather than specifying the associated physical properties
\citep{pride_relationships_2005, hubbard_hydrogeophysics_2011}. Seismic
reflection can provide high-resolution images of the CZ architecture
\citep{haeni_application_1986, bradford_depth_2002,
kaiser_ultrahigh-resolution_2009}, but its application in shallow GW is
difficult. Gradual weathering often leads to subtle or absent impedance
contrasts, hindering strong reflections. Besides, the near-surface
exhibits significant heterogeneity, noise, and the presence of surface
wave hiding reflection, necessitating dense acquisition and advanced
processing for effective imaging. SRT has been increasingly employed in
CZ exploration in the last decade \citep{olona_weathering_2010,
befus_seismic_2011, fabien-ouellet_using_2014,
holbrook_geophysical_2014, st_clair_geophysical_2015, pasquet_2d_2015,
flinchum_critical_2018, callahan_subsurface_2020, cosans_anatomy_2024},
using the inferred pressure-wave ($P$) velocity ($V_P$) as a proxy to
delineate subsurface layers (e.g.\ soil, regolith, bedrock). Borehole
data are used with these images to assign hydrofacies characteristics
\citep{befus_seismic_2011, cassidy_combining_2014, holbrook_links_2019,
guerin_borehole_2005, olona_weathering_2010, paillet_integrating_1995,
lesparre_impacts_2024}.\looseness=1

Nowadays, it is widely acknowledged that hydrofacies parameter values
(porosity, fluid saturation, and permeability) can influence seismic
properties \citep{pride_relationships_2005}. The major theory which
links hydrofacies and seismic parameters is poroelasticity proposed by
\citet{biot_theory_1956}, \citet{biot_theory_1956-2}, and
\citet{biot_mechanics_1962} as
described in the Appendix~\ref{app:poroelastic_eq}. As the S-wave
propagates primarily within the solid frame, the behaviors of P-waves
and S-waves are partially decoupled in the presence of fluids
\citep{bachrach_seismic_2000, bachrach_portable_2004}. This decoupling
implies that P-wave velocity ($V_P$) is sensitive to changes in fluid
properties, while S-wave velocity ($V_S$) primarily reflects the
elastic properties of the solid matrix. Studying both $V_P$ and $V_S$
offers valuable subsurface information, enabling more accurate
hydrofacies characterization, including porosity and fluid content in
the VZ \citep{pasquet_detecting_2015, pasquet_2d_2015} and the
piezometric surface \citep{flinchum_identifying_2020,
dangeard_river_2021}. 

While both $V_P$ and $V_S$ can be jointly estimated through SRT
\citep{turesson_comparison_2007, grelle_seismic_2009, pasquet_2d_2015},
obtaining reliable $V_S$ estimates from SRT often presents critical
challenges. These include the need for specific sources and horizontal
geophone arrays, which increases field effort, and the difficulty in
accurately identifying S-wave first arrivals during data processing. To
overcome these challenges, $V_S$ characteristics are frequently
inferred from surface-wave dispersion inversion
\citep{socco_surface-wave_2004, socco_geophysical_2010} (referred to as
Multichannel Analysis of Surface Waves, MASW). Recent studies have
extensively developed integrated approaches combining SRT and MASW to
simultaneously estimate 2D $V_P$ and $V_S$ sections along coincident
profiles from single datasets \citep{konstantaki_determining_2013,
flinchum_identifying_2020, pasquet_detecting_2015, pasquet_2d_2015,
pasquet_geophysical_2016}. With this approach, seismic methods have
become a promising imaging tool for both hydrofacies heterogeneities
and water content contrasts through the $V_P/V_S$ or Poisson's ratio
($\upsilon$) (Figure~\ref{fig:seismic_method}d). This approach,
relying on the open-source processing tools SWIP
\citep{pasquet_swip_2017, pasquet_multiwindow_2021} and PyGIMLi
\citep{rucker_pygimli_2017}, appears successful in various
hydrogeological contexts and application scales: from the laboratory on
partially saturated glass beads \citep{pasquet_small-scale_2016}, to
the field and the catchment scale \citep{pasquet_catchment-scale_2022}.
This approach has also proved to be a very promising tool for imaging
the VZ water content and the water flow paths within the river corridor
\citep{dangeard_river_2021}. In addition, it has been recently shown
that simple neural networks can effectively extrapolate 3D water table
maps from passive seismic data, using a limited number of piezometric
spatial observations \citep{cunha_teixeira_physics-guided_2025}.
Reliance on direct piezometric measured data constrains the method's
broader applicability. To address this limitation, a Transformer-based
language model inspired by neural machine translation and speech
recognition architectures has been proposed
\citep{cunha_teixeira_neural_2025}. Designed for deterministic
petrophysical inversion, it translates surface-wave dispersion data
into detailed textual descriptions of subsurface properties, including
soil composition, geomechanical parameters, and hydrogeological
conditions (basically hydrofacies). This method limits the need for
piezometric data and processes seismic data 2000 times faster than
conventional stochastic inversion methods, describing both water levels
and medium heterogeneities.

\begin{figure*}
\includegraphics{fig05}
\vspace*{2pt}
\caption{\label{fig:seismic_method}Schematic illustration from
\citet{solazzi_surfacewave_2021} and \citet{dangeard_river_2021} of (a)
a saturation-depth profile calculated with the
\citet{van_genuchten_closed-form_1980} relationship, (b) associated
porous media density $\rho_{\mathrm{pm}}(z,Sw)$, P-wave velocity $V_{p}(z,Sw)$, 
and S-wave velocity $V_{s}(z,Sw)$ profiles, and (c) body ($P$ and $S$)
and surface (PSV) waves propagating in the partially saturated soil.
(d) Sections of Poisson's ratio ($\upsilon$) estimated from seismic
data for each time step and estimated piezometric surface (colored
lines). The highest values of $\upsilon$ indicate a saturated medium,
while low values correspond to a partially or unsaturated medium The
white shaded area indicates the unconstrained area below the depth of
investigation.}
\vspace*{2pt}
\end{figure*}

A primary drawback of SRT and MASW is the manual effort required for
first-arrival picking and dispersion curve extraction. Furthermore,
quantifying the uncertainties associated with these picks remains an
active area of research \citep{dangeard_estimating_2018}. Errors are
controlled by several factors like instrument responses, poor coupling
between sources and geophones, tilt angles, near-field offset effects,
waveform changes, signal attenuation, and ambient noise sensitivity.
Data may not be reproduced accurately due to differences in source
frequency or temporal shape from location to location. Automatic source
generation cannot guarantee reproducibility in orientation, frequency
content, amplitude, and first P-arrival waveforms. These issues are not
discussed in the literature.

In recent advanced works in near-surface geophysics, the time-picking
process is generally done randomly many times (from 15 to 30 times) by
hand for each seismogram, a minimum of 15~times being necessary
\citep{dangeard_estimating_2018, dangeard_river_2021}. A standard
deviation is then estimated to quantify the picking uncertainty. Recent
deep-learning methods offer more competitive accuracy and speed in time
arrival labels compared to classical methods. Transfer and hybrid
semi-supervised deep-learning picking approaches can achieve similar
accuracies over ten times faster, especially in active seismics
\citep{huynh_near-surface_2023}. This approach is based on PhaseNet, a
U-net-based deep learning algorithm initially developed for seismic
wave picking \citep{zhu_phasenet_2018}, is improved via
Machine-Learning based corrections and extended to the context of
active seismic data.

Relative errors of approximately 5 up to 7\% are generally obtained for
P-wave first arrival for different setup configurations that are very
similar at near-surface scales \citep{pasquet_detecting_2015,
pasquet_2d_2015, bergamo_p-_2016} or at laboratory scales
\citep{bodet_elasticity_2010, bodet_small-scale_2014,
pasquet_small-scale_2016}. Besides, in those studies, dispersion curves
are also extracted from dispersion diagrams using MASW approaches.
Fundamental modes such as Rayleigh wave modes can be detected with
relative errors lower than 7--10\% by using classic random picking or
semi-supervised deep-learning approaches \citep{dai_deep_2021}.

\subsection{Petrophysical relationships for linking geophysical
observables to hydrofacies} \label{sec:petrophy}

Subsurface lithology inherently governs both hydrogeological and
geophysical data. To bridge these fields, petrophysical laws provide
quantitative links between hydrogeological properties (e.g.\ porosity,
hydraulic conductivity, water saturation) and geophysical parameters,
including seismic velocities and electrical resistivity
\citep{pride_relationships_2005, johnson_data-domain_2009,
chen_estimating_2001}. These relationships are frequently derived
through empirical calibrations based on site-specific \mbox{observations}
\citep{huisman_measuring_2003, roth_calibration_1990,
west_dependence_2001}.

These models can be refined to incorporate additional variables, such
as mineral composition, density, and various geophysical parameters.
Petrophysical models often exhibit significant variability across
different lithologies as they are highly sensitive to fluctuations in
mineral composition and pore structure \citep{mavko_rock_2009}.
Semi-empirical models, which partially integrate the geometrical and
physical characteristics of the porous medium, are widely used in
practice. The main relationships and petrophysical models used in
hydrogeophysics include:\looseness=1
\begin{itemize} 
\item Archie's law: this empirical relationship relates the electrical
conductivity of a brine-saturated rock to its porosity and water
saturation \citep{archie_electrical_1942}. Numerous studies have
adapted this law to their geological context as it is synthesized in
\citet{glover_new_2017}, \citet{glover_geophysical_2015}, and
\citet{zinszner_geoscientists_2007}. For example,
\citet{waxman_electrical_1968} and \citet{johnson_new_1986} have
extended this law to account for surface conduction in the presence of
clay, which is particularly important in low-salinity environments.
\item Gassmann's model predicts the impact of fluid saturation on the
elastic properties of rocks. It is particularly useful for the seismic
approaches \citep{gassmann_uber_1951}. 
\item Various petrophysical models and empirical relationships are used
in seismic modeling \citep{schmitt_geophysical_2015}. For instance,
Gardner's equation is an empirical relationship linking rock density to
seismic velocity \citep{gardner_formation_1974} at the near surface
scale. Hertz--Mindlin grain contact theory provides a framework for
predicting seismic velocities ($V_P$ and $V_S$) based on grain
interactions, considering factors like porosity, mineralogy, and fluid
saturation \citep{mindlin_compliance_1949}. This model was used to the
VZ by \citet{solazzi_surfacewave_2021}. The Kuster--Toks\"{o}z model
predicts the elastic properties of randomly oriented pores in rocks,
particularly in establishing relationships between lithofacies
parameters with seismic velocities \citep{mavko_rock_2009,
kuster_velocity_1974}. Additionally, the relationship of
\citet{christensen_83_2003} correlates porous media density to $V_P$
and $V_S$.
\item The empirical Kozeny--Carman relationship relates hydraulic
conductivity to porosity and specific surface area
\citep{wyllie_application_1950, carman_fluid_1937, kozeny_uber_1927}. 
\end{itemize}

These petrophysical models have been empirically adapted for CZ
applications to characterize soil water content
\citep{pasquet_catchment-scale_2022, gase_estimation_2018,
merritt_measurement_2016, pasquet_geophysical_2016, tran_joint_2014},
porosity distribution \citep{mezquita_gonzalez_quantification_2021,
callahan_subsurface_2020, flinchum_estimating_2018,
holbrook_geophysical_2014, mount_characterization_2014}, hydraulic
conductivity \citep{di_maio_combined_2015}, and fluid content
\citep{yuan_prediction_2023, holmes_application_2022}. Petrophysical
models can be derived from both laboratory measurements (e.g., core and
plug analyses, well logs) and geophysical inversions. Laboratory data
provide direct, high-resolution measurements at the core scale
\citep{schon2015physical}, {although the samples may not fully
represent in-situ conditions because their properties can be altered
during coring, extraction, and transport}. Geophysical inversions yield
indirect, spatially regularized estimates that integrate over larger
volumes and depend on prior assumptions and inversion constraints
\citep{tarantola_inverse_2005,aster_parameter_2018}. Both sources are
widely used and typically complement each other in subsurface
characterization \citep{pride_relationships_2005,linde_joint_2016}.\looseness=1

A key challenge arises when applying petrophysical laws to these
geophysical inverted fields. Geophysical inversions produce spatially
correlated models due to both the physics of data acquisition and the
regularization applied during inversion. Typically, these fields
exhibit long-range continuity in the horizontal direction and
short-range continuity in the vertical direction, reflecting geological
layering and the geometry of surface measurements
\citep{chiles_geostatistics_2012,
isaaks_introduction_1989,grana_geostatistical_2022}. In addition,
poorly constrained petrophysical relationships themselves remain a
major source of uncertainties in models
\citep[e.g.][]{seo_improved_2023, tso_field_2019, brunetti_impact_2018}.\looseness=1

\subsection{Geostatistical methods for subsurface characterization}

Integrating diverse data types requires the application of robust
geostatistical algorithms. While high-resolution spatial studies often
rely on datasets with limited 2D extent, there is a growing need to
incorporate 3D or 4D (timelapse) dimensions. To address these
complexities, several powerful geostatistical frameworks are available.
The package gstlearn \citep{renard_gstlearngstlearn_2025}, developed by
Mines Paris, offers a comprehensive suite of tools for such
multivariate challenges, complementing other established packages like
geoR \citep{ribeirojr_geoR_2025}, gstat \citep{pebesma_gstat_1998}, or
GEONE (DeeSse) \citep{zhexenbayeva_deesse_2024}. A key challenge is to
construct a three-dimensional model providing the hydrofacies
distribution by assigning probabilities to the geophysical property
values associated with a particular medium
\citep{gottschalk_integrating_2017, zhu_improved_2016}. 

Hard data, or quantitative data, refers to information that is
measurable and verifiable. It is a type of data that is measured and
can be analyzed statistically. In contrast, soft data and qualitative
data are mostly derived indirectly (for example, related to the
variable of interest), are usually more numerous and less accurate than
hard data because of their indirect nature. One, two, or
three-dimensional spatial information derived from geophysical imaging
and hydrogeological interpretation (e.g.\ cross-sections or calibrated
flow models) are inherently uncertain and therefore considered as soft
data. Combined with geological logs (hard data), they could be used to
perform conditional geostatistical simulations.\looseness=1

Ordinary kriging has also been applied in critical zone studies to
interpolate CZ architecture, as demonstrated by a 3-D seismic velocity
model constructed beneath a soil-mantled granite ridge in the Laramie
Range, Wyoming, using 25 seismic refraction transects
\citep{flinchum_critical_2018}. This geostatistical approach enabled
the delineation of subsurface layers and improved understanding of deep
CZ structure. Similar kriging-based strategies have been implemented to
map critical zone compartment geometry, such as in the Strengbach
watershed (Vosges Mountains, France) \citep{lesparre_impacts_2024}, and
in the Quiock watershed (Guadeloupe archipelago, France)
\citep{pasquet_catchment-scale_2022}.

The multiple-point simulation (MPS) approach developed by
\citet{strebelle_conditional_2002} can be used to provide simulated
facies maps. In practice, it directly uses empirical distributions
inferred from one or several training images that can be derived either
from similar outcrop observations, expert knowledge, or geophysical
data \citep{caers_stochastic_2003, linde_geological_2015}. Among the
MPS methods, Direct Sampling is often employed since it is very
flexible and computationally efficient
\citep{meerschman_practical_2013}. It has been recently extended for
the simulation of multi-resolution patterns
\citep{straubhaar_multiple-point_2020}, or for handling complex
conditioning data such as inequality constraints
\citep{straubhaar_conditioning_2021}. These methods reproduce complex
geometries such as channels, meanders, or lenses while preserving the
relations between the facies.

An alternative is the Plurigaussian Geostatistical Simulation (PGS)
approach. PGS has long been acknowledged as a practical tool for
assessing the impact of subsurface heterogeneity on field and basin
scale flow and transport processes \citep{armstrong_plurigaussian_2011,
le_loch_improvement_1994, abzalov_plurigaussian_2023}. The outcome of
MPS and PGS simulation methods is a set of 3-D equiprobable
realizations of hydrofacies that are conditioned by the hard data and
can be used to derive the probability of occurrence for each
hydrofacies at a given location.

Process-based or hybrid stochastic/genetic models, such as Flumy, can
be employed to generate geologically consistent prior realizations.
These models more faithfully capture sedimentary dynamics while
remaining strongly conditionable to geophysical and hydrological data.
In this way, Flumy-type priors can complement geostatistical
interpolation by providing more realistic subsurface structures that
serve as a robust basis for hydrogeophysical inversion across different
dimensionalities (1D--3D). Moreover, \citet{hermans_uncertainty_2015}
demonstrated how training-image-based priors can be updated or
falsified using electrical resistivity tomography (ERT) data, improving
the consistency between geological scenarios and geophysical
observations. Building on this, \citet{neven_stochastic_2022}
introduced a multi-fidelity inversion framework that combines
low-fidelity geologically consistent stochastic models with
high-fidelity GW flow simulations, enabling efficient joint inversion
that preserves geological realism while reducing computational cost.

While such methods are suitable for the simulation of heterogeneous
media, the addition of geophysical and hydrological data (secondary
information) makes them potentially more powerful. Bayesian sequential
simulation method (BSS) from \citet{doyen_bayesian_1996} provides an
attractive alternative for simulating jointly primary and secondary
variables, without requiring any linear or specific parametric
relationship between the variables. The method is based on PGS, with
the addition that the conditional distribution estimated by kriging is
updated by the joint distribution of the primary and secondary
variables in a Bayesian framework. While BSS is recognized as a
powerful and flexible geostatistical tool
\citep{chen_estimating_2001,chiles_fifty_2018}, it also provides
significant improvements. In particular, the simultaneous reproduction
of: (i) the variance and the fine-scale structure, and (ii) the
relationship between the primary and the secondary variables. It is
particularly adapted for the hydrogeophysical applications, in order to
deal with the smooth behavior of the secondary variable due to the
regularization procedure involved in the geophysical inversion
procedures.\looseness=1 

The traditional geostatistical framework assumes that the spatial
characteristics of the variables are defined and valid over the whole
field. It fails at reproducing local changes. For that sake, the new
branch of non-stationary Geostatistics has been recently launched which
offers the possibility to consider any parameter of the spatial
characteristics (e.g.\ range, anisotropy orientation, variogram sill) as
defined locally \citep{renard_gstlearngstlearn_2025}. To be fruitful,
this technique requires additional information such as a proxy driving
the anisotropy orientation (following a meander direction for example)
or a proxy for the variance (following the river flow for example)
\citep{fouedjio:tel-01139460}.\looseness=1 

\section{Geophysical and hydrological wedding}

Historically, a disciplinary divide has persisted between geophysicists
and hydrogeologists, largely driven by diverging research priorities.
While geophysicists have traditionally focused on advancing
instrumentation, inversion algorithms, and interpretation frameworks,
hydrogeologists are deeply driven by the development of predictive
models to simulate complex hydrological processes within the CZ. This
disconnect often arises from a difference in application:
hydrogeologists frequently seek to optimize their numerical models
using geophysical data as a spatial constraint, whereas geophysicists
have historically prioritized the methodological validation of the
signal itself over the broader hydrological outcomes
\citep{robinson_advancing_2008}. The emerging fields of hydrogeophysics
\citep{binley_emergence_2015} and CZ sciences aim to bridge the gap
between geophysical and hydrological data to explore their
interrelationships. Both disciplines seek to understand subsurface
processes using the same Laplace's equations. Researchers are
increasingly focused on correlating equation parameters, converting
geophysical measurements into hydrological features and vice versa, to
develop more precise, data-driven models.

\subsection{Conceptual framework for data and model integration}

Integration of geophysical and hydrological models requires robust
inversion frameworks capable of reconciling sparse and indirect
observations with inherent complexity of subsurface dynamics. The
inversion process seeks a parameter distribution that minimizes a
misfit function (i.e.\ the discrepancy between observed and synthetic
data); in order to identify the physical properties that govern
processes such as GW flow or electrical resistivity. In practice, these
inversion strategies are generally categorized into deterministic and
stochastic approaches.\looseness=1 

As illustrated in Figure~\ref{fig:workflow_inv_pb}, the methodologies
for integrating hydrogeophysical data follow two main workflows: (i)
sequential hydrogeophysical inversion (SHI), where data are processed
through separate and/or successive stages; (ii) processes based
inversion (PBI), which solves for the underlying transient
process-based parameters directly \citep{hellman_structurally_2017,
herckenrath_sequential_2013, hinnell_improved_2010,
wagner_overview_2021,wagner_overview_2021,Linde2017AdvancesHydrogeophysics,
hubbard_hydrogeophysics_2011}.

\begin{figure*}
\includegraphics{fig06}
\vspace*{2pt}
\caption{\label{fig:workflow_inv_pb}(a) Sequential hydrogeophysical
inversion (SHI): geophysical data are independently inverted and
interpreted to inform the structural framework of the hydrogeological
model, which is then calibrated using hydrological observations. (b)
Process-based inversion (PBI): hydrological and geophysical forward
models are dynamically linked. Multiple geophysical and/or hydrological
datasets are integrated within a unified inversion framework, with
coupling terms or petrophysical constraints. Arrows represent the flow
of information as well as iterative parameter and condition updating,
continuing until a minimum data and/or constraint misfit is achieved.}
\vspace*{2pt}
\end{figure*}

\subsubsection{Deterministic inversion: local optimization}

Deterministic inversion methods are based on the computation of the
gradient of the misfit function using local iterative linearization
approaches. In these common local approaches, to reduce the space of
possible solutions, a regularization procedure is applied
\citep{menke_describing_2012} often imposing smoothness or closeness to
a prior model. If the initial solution is too far from the global
solution, the inversion can be stopped in a local minimum or in an area
where the cost function does not significantly change. Linearized
approximations like those used by \citet{cockett_efficient_2018} remain
computationally efficient, but may fail in highly nonlinear or
ill-posed settings \citep{dagasan_using_2020}. Moreover, local
approaches are not suitable for estimating reliable uncertainties.

\subsubsection{Stochastic inversion and uncertainty quantification}

In contrast, stochastic inversion approaches aim to explicitly quantify
uncertainty by sampling the posterior distribution of model parameters
\citep[e.g.][]{vrugt_hydrologic_2013, rings_coupled_2010,
hinnell_improved_2010, chen_coupled_2003,
sambridge_monte_2002,hubbard_hydrogeological_2000}.
\citet{dmowska_chapter_2014} give a comprehensive description of
existing methods for estimating uncertainties and illustrate that a
posteriori probabilistic distribution sampling of model parameters
within a Bayesian framework is much more suitable. As a matter of fact,
importance sampling enables the cost function to be sampled
proportionally to the probability \citep{tarantola_inverse_2005}. This
sampling is generally carried out using Markov chain Monte Carlo
(MCMC)-based algorithms
\citep{mosegaard_monte_1995,sambridge_monte_2002}. State-of-the-art
inversion \citep[e.g.][]{aster_parameter_2018,zhou_inverse_2014,
tarantola_inverse_2005} and data assimilation (e.g., ensemble smoother,
ensemble Kalman filers and its variants)
\citep[e.g.][]{mcaliley_application_2022,jiang_dart-pflotran_2021,
chen_application_2013} are generally used to ingest the hydrological
and geophysical data.

\subsubsection{Misfits definition}

The combination of various misfits related to dynamic data, such as
water level or temperature, with those of spatial data, including
travel time, water level geometry, or electrical resistivity, is not
straightforward. In the sequential hydrogeophysical inversion, the
objective functions are estimated independently, whereas in the joint,
coupled or process-based inversion, the objective function can assume
various forms \citep{mboh_coupled_2012, hinnell_improved_2010,
rubin_hydrogeophysics_2005} (see section below). The misfit of each
measurement can be weighted in the objective function according to the
measurement uncertainty and sensitivity of each method.

\subsection{Static hydrofacies characterization by Sequential
hydrogeophysical Inversion approaches (SHI)}

In a SHI, geophysical data is separately inverted to estimate the
distribution of a single geophysical property (e.g., a map of
resistivity). Estimated geophysical property maps are subsequently
employed to delineate hydrofacies structures, including the number,
geometry, and spatial distribution of facies boundaries of the
hydrogeological FWD (Figure~\ref{fig:workflow_inv_pb}a) or state
variable (e.g., saturation or salinity from resistivity)
\citep{binley_emergence_2015, singha_advances_2015,
razafindratsima_hydrogeological_2014, bendjoudi_riparian_2002,
hubbard_hydrogeological_2000, rubin_mapping_1992}. In these workflows,
geophysical and hydrological models are processed sequentially rather
than iteratively. The objective functions of the geophysical and
hydrological models are not linked; each model is calibrated
independently without direct coupling between their respective misfit
functions. The geophysical data are first interpreted independently to
infer subsurface structure, which then informs the hydrogeological FWD
(Figure~\ref{fig:workflow_inv_pb}a). After defining the inversion
framework, GW flow and transport simulations, whether stochastic or
deterministic, are conditioned on prior estimates of hydrofacies
parameters such as hydraulic conductivity or porosity. These prior
ranges, which help constrain parameter uncertainty, are typically
derived from laboratory analyses, pumping tests, or empirical
correlations with well-log or geophysical data
\citep{hyndman_coupled_1994, ahmed_combined_1988}. In addition to
defining prior ranges, certain measurements, such as pumping tests, can
also be incorporated directly into the inversion process. The
calibration of the GW model parameter is then carried out while keeping
the conceptual model fixed \citep{kollet_integrated_2006}. The earliest
study that used this approach was conducted by
\citet{ahmed_combined_1988}, which used the co-kriging of measured
transmissivity, specific capacity, and electrical resistivity to
elaborate transmissivity maps.

Geophysical joint inversion involves integrating multiple geophysical
data types to improve the characterization of subsurface properties.
This approach is particularly beneficial in complex heterogeneous
environments \citep{gallardo_joint_2007}. To reduce such drawbacks,
joint inversion methods combining SRT and ERT simultaneously have been
developed at the near-surface scale \citep{doetsch_imaging_2012,
gallardo_joint_2007, gallardo_joint_2004} and could be extended to
hydrogeophysics problems. Different constraints linking directly the
$P$ and/or $S$ velocities to resistivity, such as petrophysical models,
cross-gradients (similarity constraints) can be introduced in the
misfit function to reduce the models' uncertainties during the
inversion process. Furthermore, a priori geological information such as
lithological boundaries could also be inverted during the inversion
process at the same time as the property values or separately by
introducing level-set constraints as in \citet{giraud_structural_2021}, and
\citet{zheglova_multiple_2018}. Recent studies have adopted multiple-point
geostatistics or geological scenario modeling to generate ensembles of
plausible hydrofacies configurations
\citep{hermans_uncertainty_2015,neven_novel_2023}. These scenarios are
subsequently evaluated by comparing simulated responses with observed
geophysical data, often through a probabilistic falsification
procedure, providing a more robust basis for structural inference. By
extending those approaches to hydrogeophysics, more particularly to
seismic data and ERT, it could be possible to reconstruct geological
features compatible with different physics but also to retrieve
physical properties that are not necessarily correlated. If $V_P$
seismic models using ray-tracing-based inversions and $V_S$ seismic
models using wave dispersion curves inversions can be obtained, they
can also be improved by using full waveform inversion at the
hydrogeophysical scale \citep{groos_application_2017,
groos_challenges_2014,liu_critical_2022}. {For example,
\citet{eppinger_2d_2024} developed a full-waveform tomography workflow
to investigate 2D weathering patterns in the critical zone, applied in
the Medicine Bow National Forest (USA). Their results are promising,
but they also point out that extending the approach to 3D would require
substantially greater computational resources.}\looseness=1

Besides, joint inversions can be performed using various combinations
of physics, such as different heat tracing, ERT, Ground Penetrating
Radar, SRT, and MASW, and constrained by petrophysical models relating
parameters that characterize the medium under study. Those
petrophysical constraints allow not only to reduce the number of free
parameters but also to obtain a final inverted model compatible with
the different physics \citep{qin_indirect_2024,
steiner_quantitative_2022, wagner_overview_2021,
wagner_quantitative_2019}.

As long as the inversion is deterministic, a key limitation of
geophysical joint inversion approaches is their reliance on multiple
geophysical datasets to define a single groundwater conceptual model,
which may not fully capture subsurface heterogeneity. In addition, this
approach can create discrepancies related to state variables such as
water saturation or pressure, and it does not account for the transient
behavior of these variables, which is critical in dynamic groundwater
systems. If the geometry was fixed before inverting for hydrogeological
parameters, and not correctly estimated, the hydrogeological parameters
will likely have to be incorrectly identified during inversion in order
to compensate for these initial errors \citep{de_marsily_current_1998}.
The outputs of GW simulations, such as saturation, hydraulic heads,
fluxes, or transport behavior, do not influence the geophysical
inversion, while the hydrofacies parameters influence geophysical data.
This issue becomes particularly critical under transient conditions,
where the dynamic evolution of water and heat fluxes provides valuable
information that could help refine both structural and parametric
interpretations of the geophysical data. Recent studies have adopted
multiple-point geostatistics or geological scenario modeling to
generate ensembles of plausible hydrofacies configurations
\citep{hermans_uncertainty_2015, enemark_hydrogeological_2019,
enemark_incorporating_2024}. 

\subsection{Toward transient hydrogeophysical inversions}
\label{sec:transient}

The CZ requires time-lapse geophysical surveys to move beyond static
imaging and effectively monitor and quantify dynamic subsurface
processes over time, thus transitioning geophysics from a field
primarily concerned with static structural imaging to one adept at
monitoring and quantifying dynamic physical and biogeochemical
processes. This continuous monitoring has demonstrated efficacy in
identifying transient processes in the CZ, such as water infiltration
\citep{daily_electrical_1992, wehrer_characterization_2016,
hu_advanced_2023, blazevic_time-lapse_2020}, variations in water
storage \citep{galibert_quantitative_2016}, flowpaths identification
\citep{nyquist_testing_2018} GW recharge, and the propagation of
thermal fronts \citep{shariatinik_ert_2024, lesparre_4d_2019,
hermans_quantitative_2015}. The three first methods are based on the
SHI approaches.

\subsubsection{Absolute inversion}

Absolute inversion is the method commonly used to interpret time-lapse
datasets. Every geophysical dataset is inverted independently.
Transient subsurface changes are then identified by comparing the
resulting models for each time step, either on the basis of absolute
values or relative differences \citep{herckenrath_sequential_2013,
LaBrecque_1995,daily_electrical_1992, slater_cross-hole_2000}. This 
method is computationally efficient and simple. However, some studies
have shown that this approach can result in unrealistic geophysical
parameters changes because it is sensitive to data noise and does not
directly account for time correlations or hydrological constraints
\citep{miller_application_2008, loke_use_2022, johnson_time-lapse_2015,
labrecque_difference_2001}.  To address the limitations of absolute
inversion, several hydrogeophysical strategies have been developed, as
discussed in the following subsections.

\subsubsection{Difference inversion}

Differencing and ratio methods focus directly on changes between
datasets, rather than reconstructing each model separately, allow to
decrease these effects \citep{labrecque_difference_2001}. The
difference inversion enhances sensitivity to temporal changes by
inverting differences in data $ \delta d_t = d_t - d_0 $, and model
parameters $\delta m_t = m_t - m_0$, where $m_0$ is the baseline model.
The inversion is typically linearized around the baseline using a
Jacobian matrix $J_{ij} = ({\partial d_{t,\mathrm{sim},i}})/({\partial
m_{t,j}})$ which quantifies the sensitivity. Introduced by
\citet{labrecque_ert_1996}, it can also be computationally efficient if
the Jacobian is reused across time steps
\citep{labrecque_difference_2001}. This approach has proven valuable in
resolving subsurface changes with time-lapse ERT during heat-tracing
experiments and solute transport, as shown by
\citet{hermans_covariance-constrained_2016} and
\citet{bretaudeau_time-lapse_2021}, who applied covariance-based
constraint of model differences and double-difference inversion schemes
to improve sensitivity and reduce artifacts.\looseness=1

\subsubsection{Time-constrained inversion} 

The time-constrained inversion introduces temporal smoothness
constraints to stabilize the inversion across time steps. It minimizes
a cost function that balances two terms: (i) the data misfit, measuring
agreement with the observed responses, and (ii) a temporal
regularization term, penalizing abrupt parameter changes over time. The
optimization is frequently solved using Gauss--Newton iterations, with
the Jacobian matrix $J_{ij} = ({\partial d_{t,\mathrm{sim},i}})/({\partial
m_{t,j}})$ quantifying the sensitivity of simulated data to parameter
updates. 

Several strategies have emerged. Sequential approaches invert each time
step independently, using a baseline reference model (typically from
the initial survey, $t_0$) to compute resistivity changes at subsequent
time steps \citep{miller_application_2008}. The baseline model ($t_0$)
is typically acquired under reference conditions (e.g., before
hydrological perturbations) and serves as the structural and
resistivity reference for all subsequent surveys. In this framework,
each survey is inverted separately, and time-lapse changes are analyzed
by comparing the inverted models post-inversion. While computationally
efficient, this approach can introduce noise amplification and
inconsistencies between surveys, as each inversion is subject to
different regularization choices and data quality variations, which can
lead to anomalous artifacts \citep{lesparre_4d_2019}.

In contrast, joint inversion frameworks have been proposed, in which
all time steps are inverted simultaneously while imposing temporal
smoothness constraints \citep{kim_four-dimensional_2013,
karaoulis_4d_2014, wilkinson_windowed_2022, zhou_joint_2025,
karaoulis_time-lapse_2013, doetsch_structural_2010}. 
These methods improve model stability and reduce artifacts in dynamic
environments like riverbeds \citep{mclachlan_limitations_2020}, during
solute transport \citep{karaoulis_4d_2014}. However, excessive
smoothing may obscure real changes, and coarse inversion meshes may
limit resolution \citep{liu_rapid_2020, hermans_quantitative_2015}.

One promising strategy is to embed geostatistical information that
characterizes the spatial correlation structure of the subsurface. For
example, in an alluvial context, \citet{jordi_geostatistical_2018}
constructs regularization operators based on an exponential covariance
model that describe the heterogeneities of the subsurface. Applications
of this approach to layered alluvial systems are demonstrated by
\citet{palacios_time-lapse_2020} and
\citet{arboleda-zapata_time-lapse_2025}.

\subsubsection{Transient process-based inversion: joint inversion or
fully coupled hydrogeophysical inversion}

Process-based hydrogeophysical inversion refers to approaches that
dynamically couple hydrological and geophysical FWDs, enabling
simultaneous data integration and bidirectional interaction during the
inversion process. In time-lapse applications, these methods
incorporate temporal evolution of subsurface parameters and state
variables \citep{hermans_advancing_2023}. The hydrogeophysical
literature uses various terms, including ``coupled hydrogeophysical
inversion'', ``joint hydrogeophysical inversion'', ``process-based
inversion'' \citep{wagner_overview_2021}, and ``fully-coupled
inversion'' \citep{hubbard_hydrogeophysics_2011}, often interchangeably
to describe approaches that integrate hydrological and geophysical data
\citep{linde_joint_2016, finsterle_joint_2008}. While terminology
varies across the community, these approaches share the common goal of
coupling hydrological and geophysical models through petrophysical
relationships to avoid the biases inherent in sequential approaches
\citep{daylewis_applying_2005, linde_joint_2016,
herckenrath_sequential_2013}. The key differences lie in implementation
strategies---specifically, how petrophysical relationships are
incorporated into the inversion framework---rather than in formal
nomenclature.

In one common implementation (often termed coupled hydrogeophysical
inversion), the hydrological FWD provides state variables such as water
saturation, salinity, or temperature, which are converted into
geophysical parameter distributions through petrophysical models
(Figure~\ref{fig:workflow_inv_pb}b)
\citep{hubbard_hydrogeophysics_2011, huisman_hydraulic_2010,
linde_joint_2016, wagner_overview_2021, bascur_hybrid_2025}. The
parameters of the hydrological FWD are updated by minimizing a combined
objective function that accounts for both hydrological and geophysical
misfits \citep{mboh_coupled_2012, hinnell_improved_2010}. A key
advantage is that geophysical inversion is guided by hydrological
process models, with petrophysical-hydrological regularization avoiding
the potential biases arising from traditional geophysical
regularization, such as smoothness and smallness constraints
\citep{daylewis_applying_2005, linde_joint_2016}.

An alternative implementation strategy (often termed joint
hydrogeophysical inversion) processes multiple datasets simultaneously,
enforcing structural or petrophysical consistency across geophysical
and hydrological FWDs (Figure~\ref{fig:workflow_inv_pb}c)
\citep{herckenrath_sequential_2013, steklova_joint_2017}. This approach
employs a modular framework that integrates the parameters of
hydrological and geophysical models into a unified vector. The
inversion process minimizes an integrated objective function comprising
three weighted components: geophysical data misfit terms ($\Phi_{SWD}$,
$\Phi_{TT}$, $\Phi_{\rho}$), hydrogeological data misfit terms
($\Phi_{pw}$, $\Phi_{Sw}$, $\Phi_{T}$), and an explicit coupling
constraint misfit ($\Phi_{\mathrm{coup}}$). The most common coupling
mechanisms are petrophysical relationships
\citep{wagner_overview_2021,steklova_joint_2017}. This strategy differs
from the previous approach in that both hydrological and geophysical
parameters are simultaneously estimated, with petrophysical
relationships serving as explicit constraints rather than forward
transformations.

A critical challenge in both approaches is the appropriate weighting of
different misfit terms to balance the influence of diverse data types
with varying units, sensitivities, and uncertainties
\citep{linde_joint_2016}. Common approaches include equal weighting of
normalized misfits, uncertainty-based weighting (where weights are
inversely proportional to data variance or covariance matrices), and
data error-based weighting \citep{mboh_coupled_2012,
kowalsky_joint_2006, linde_joint_2016}. \citet{mboh_coupled_2012}
demonstrated that uncertainty-based weighting using the inverse of data
covariance matrices provides a statistically rigorous framework for
balancing electrical resistance and hydrological measurements in
coupled inversion. \citet{christensen_testing_2016} showed that
propagating geophysical data uncertainties into hydrological model
calibration significantly affects prediction error, highlighting the
importance of proper error-based weighting strategies.  More
sophisticated strategies have emerged, including: (i) adaptive
weighting schemes that adjust weights during the inversion process
\citep{steklova_joint_2017}, (ii) ensemble-based data assimilation
approaches where weighting is handled through estimated data-model
covariances rather than fixed scalar weights \citep{neven_novel_2023},
(iii)~multi-objective Pareto optimization that avoids pre-selecting
scalar weights by producing trade-off solution fronts
\citep{danek_pareto_2019}, (iv) Bayesian frameworks where weights
naturally emerge from data and prior covariance matrices
\citep{dettmer_probabilistic_2024}. However, weight selection remains
partly subjective and can significantly influence inversion results,
particularly when data quality or petrophysical uncertainty varies
across datasets \citep{commer_hydrogeophysical_2013,
finsterle_itough2_2017}.

In alluvial aquifer settings, \citet{herckenrath_sequential_2013}
demonstrated that this approach improves parameter estimation and model
consistency compared to sequential hydrogeophysical inversion when
high-quality petrophysical constraints are available, though at a
higher computational cost. Recent applications highlight both the
opportunities and challenges of hydrogeophysical inversion. For
example, \citet{pleasants_hydrogeophysical_2022},
\citet{pleasants_comparison_2023} and \citet{boyd_coupled_2024} used
time-lapse ERT data in mountainous regions to improve groundwater flow
models. In these studies, simulated water saturation was first
converted to resistivity using petrophysical relationships (e.g.,
Archie's Law). A weighted misfit function was then employed to compare
these modeled resistivities with ERT data and to balance the influence
of the different data types during calibration.
\citet{bascur_hybrid_2025} introduced a Hybrid Bayesian Inversion
framework for CHI using single-time 2-D ERT data from an unconfined
alluvial aquifer. Their approach jointly integrates saturated GW FWD
modeling with geophysical observations via uncalibrated petrophysical
models, with weighting naturally emerging from the Bayesian framework
through data and prior covariance matrices. Tested on both synthetic
and experimental datasets, this framework provides efficient
uncertainty quantification and is particularly suited to estimating
water fluxes in heterogeneous aquifers. To address computational
challenges, \citet{neven_stochastic_2022} proposed a stochastic
multi-fidelity framework integrating Ensemble Smoother with Multiple
Data Assimilation (ES-MDA) and multi-fidelity modeling. Applied to both
synthetic and real datasets, this approach improved accuracy while
significantly reducing computational costs.

The recent development of a fully coupled GW--ERT FWD using pore water
pressure as the coupling variable \citep{salas-ariza_permeability_2024}
represents a promising methodological advancement. A sensitivity
analysis of this model was performed, though full inversion using the
coupled framework in real field applications had not yet been
demonstrated at the time of publication.

Both implementation strategies are typically applied to sites with
relatively simple geometries, low heterogeneity, and limited
multiphysical calibration datasets. Future research is needed to
address the definition of coupling terms between different FWDs, the
assignment of weights in objective functions for various data types,
and the evaluation and quantification of petrophysical uncertainty.
While coupled hydrogeophysical inversion has demonstrated success in
controlled settings and moderately complex field sites, extending these
approaches to highly heterogeneous systems with complex geometries
remains challenging. 

\section{Remaining challenges and future directions}

Shallow groundwater processes in the critical zone occur within
strongly heterogeneous environments, where heat, water, and solute
fluxes vary in space and time \citep{huggins_groundwater_2023}. Surface
interfaces such as SW--GW boundaries and the vadose zone are especially
heterogeneous (see Appendix~\ref{app:geology}), making it difficult to
resolve the associated four-dimensional water and heat fluxes. Despite
recent advances in hydrogeophysics and modeling, accurately
characterizing and predicting these dynamics remains challenging,
particularly in terms of data interpretation and process representation
\citep{binley_emergence_2015, romeroruiz_review_2018,
ledford_connecting_2022, whiteley_geophysical_2019,
dangeard_river_2021, hermans_advancing_2023}. A key difficulty arises
from the mismatch between data types: hydrological measurements
typically provide high-accuracy, high-temporal-resolution information
but only at sparse spatial locations, whereas geophysical observations,
{such as ERT or SRT}, offer spatially extensive coverage but at coarser
spatial resolution and with lower temporal resolution and greater
uncertainty. Reconciling these datasets, therefore, requires
multivariate methodologies capable of integrating information with
different spatial supports and error structures
\citep{wackernagel_multivariate_2003}.

Together, these complexities limit our ability to fully resolve
four-dimensional patterns of shallow GW flow, to quantify uncertainty,
and to design models and monitoring strategies that are both
representative and computationally tractable. Addressing these
limitations requires careful consideration of: (i)~the balance between
model complexity and representativeness; (ii) approaches for managing
parameter dimensionality and uncertainty; (iii) integration of
multi-dimensional geophysical and hydrological datasets; (iv) treatment
of transient parameters and state variables; (v) model-guided
strategies for optimal monitoring design; and (vi) scaling frameworks
that enable transition from hillslope to catchment domains while
accounting for change-of-support issues.

The following subsections discuss these remaining challenges in detail
and outline opportunities for advancing the characterization and
modeling of shallow GW systems in the CZ. Highly instrumented
observatories such as CZNet and eLTER networks (e.g., OZCAR, TERENO)
offer unique opportunities in this context. These sites are not only
geologically well characterized but also benefit from dense spatial
networks and extended time series of hydrological measurements, which
are critical for defining robust initial and boundary conditions in
numerical models \citep{gaillardet_ozcar_2018}.

\subsection{Multi-dimensional geophysical and hydrological data}

A fundamental difficulty in hydrogeophysical characterization lies in
the integration of complementary yet dimensionally and temporally
distinct data types. On the one hand, 2D ($x,y$) or ($x,z$) and 3D
geophysical surveys provide spatially extensive but temporally discrete
information. On the other hand, hydrological observations (e.g.,
temperature, water level, discharge, or high-frequency DTS
measurements) offer high temporal resolution but are limited to point
or line scales. Combining these heterogeneous datasets to reveal
spatiotemporal patterns, such as water-table dynamics,
temperature-field evolution, soil-moisture variations, and water-flux
behaviors, therefore remains methodologically challenging.

Accurately defining initial and boundary conditions, as well as
transient parameters, also remains a key challenge in GW flow and
transport models within the CZ, particularly in unsaturated,
heterogeneous environments and data-scarce regions. Transient dynamics
are frequently simplified in GW FWD due to insufficient temporal
resolution or overly idealized boundary conditions. Key hydrological
and geophysical variables are often observed at sparse and irregular
intervals across space and time. Predicting their 4D distribution,
therefore, requires robust interpolation or modeling frameworks capable
of bridging these observational gaps.

The sensitivity of parameter estimation to initial conditions has been
widely studied in hydrology, with evidence showing that different
initial parameter values can lead to significantly varied model
outcomes highlighting that initial conditions are as critical as the
calibration process itself \citep{liu_simultaneous_2009}. Recent
studies have demonstrated that geophysical data can be used to estimate
initial states and boundary conditions of GW models. For instance,
integrating seismic time-lapse data with hydrological time series
enables the reconstruction of the initial shape of the water table and
the boundary conditions, thereby facilitating the simulation of water
flow paths in river corridors \citep{dangeard_river_2021}. Another
approach, as developed by \citet{bascur_hybrid_2025}, is to {suggest}
hybrid Bayesian inversion frameworks to integrate single-time ERT data
and infer initial water table depth along with hydraulic conductivity,
even under petrophysical uncertainty. Similarly,
\citet{pleasants_hydrogeophysical_2022} derive hydraulic parameters and
initial conditions in a hillslope-scale model from time-lapse ERT,
showing that coupled hydrogeophysical inversion improves physical
consistency and model initialization.

However, these approaches still face limitations when data coverage is
sparse or when geophysical measurements are irregularly distributed in
space and time. This calls for a multivariate methodology with
different supports of information
\citep{wackernagel_multivariate_2003}. The traditional geostatistical
framework usually implies that the spatial characteristics of the
variables are defined and valid over the whole field. These constraints
have been addressed by introducing the local geostatistics concept
\citep{renard_modeling_2013,thiesen_assessing_2022} where the
parameters of the spatial characteristics (variograms) can vary over
the domain. This improvement yet implies the use of a local approach
(through a compulsory moving neighborhood).  More recently, this
flexibility has been reinforced by maintaining consistency even at
large scale, through the stochastic partial differential equations
(SPDE) \citep{renard_inference_2019, clarotto_spde_2024}. This
methodology presents two advantages: it allows handling complex
non-stationarity where any variogram parameter could follow information
provided by auxiliary information; it also enables efficient processing
of large amounts of information, ignoring any neighborhood limitation
\citep{lindgren_explicit_2011}. 

Finally, the SPDE model incorporates physics within a stochastic
framework. It is well defined when processing a variable in a space and
time domain. For example, \citet{clarotto_spde_2024} develops an
advection--diffusion version of SPDE to define a large class of
spatio-temporal models. \citet{grimaud_kriging_2025} demonstrate the
efficiency of SPDE-based kriging in modeling nonstationary anisotropies
for predicting alluvial thickness, yielding more geologically realistic
patterns and lower uncertainties than traditional stationary
approaches. SPDE could offer more accurate representations of
asymmetrical space--time processes, such as those found in
hydrogeophysical phenomena. Incorporating time into PGS is essential to
accurately represent water content variations over time and interpolate
the hydrofacies in 4D. The time variations of the hydrofacies are
linked to water saturation fluctuations. A water level mapping
geostatistical framework that accounts for SW--GW connectivity has for
instance been elaborated by \citet{maillot_technical_2019}. The method
is conditioned by `hard' data (hydraulic head) and ``soft'' data (dry
wells) as conditioning points. The next challenge is to incorporate the
2D water level produced by the geophysical methods.

\subsection{Coupled properties inversion}

{Hydrogeophysical process-based inversion must reconcile not only
contrasting spatial and temporal scales but also the fundamentally
different physical meanings of hydrological and geophysical
observables.} Hydrological observables (e.g., pressure head, solute
concentration, temperature) and geophysical observables (e.g.,
resistivity, seismic velocity) respond to distinct physical mechanisms,
such as flow and transport, electromagnetic or elastic properties,
which do not necessarily share a common characteristic length or
wavelength, nor the same way of transferring, propagating, diffusing
or, more generally, reflecting the medium under study
\citep{Linde2017AdvancesHydrogeophysics}. Establishing meaningful
correlations therefore requires identifying parameters that are jointly
sensitive to both domains, such as porosity, water saturation, or
salinity. However, these parameters represent effective quantities
averaged over distinct scales, meaning that one inferred parameter may
not describe the same physical property or structural feature across
domains. This scale discrepancy complicates the establishment of robust
correlations between hydrological and geophysical models. Coupling
should therefore be restricted to parameters or processes that share
comparable scales of sensitivity, while recognizing that others remain
uncorrelated due to fundamental differences in resolution and physical
response. Although several studies have proposed petrophysical or
constitutive laws to link hydrological and geophysical parameters,
these relationships generally fail to resolve inherent scale
discrepancies, as they are often defined at the laboratory or local
scale and lose validity under heterogeneous field conditions
\citep{binley_resistivity_2020}.

To address this, the field is moving toward physics-informed machine
learning \citep{cunha_teixeira_physics-guided_2025} and probabilistic
upscaling, or U-net based deep learning joint multiphysics inversion
models {\citep{hu_deep_2023, ren_joint_2024,oh_deep_2024}} where models
learn effective relationships from multi-domain datasets rather than
relying on rigid constitutive laws \citep{neven_stochastic_2022}, with
or without structural similarity constraints (such as cross-gradients,
level-set structural geometries,\,{\ldots}). These approaches combine
physical constraints with data-driven flexibility, enabling improved
parameter estimation and uncertainty quantification in heterogeneous
environments.

\subsection{Choosing model spatial dimensionality in coupled
hydrogeophysical inversion}

{Determining the appropriate spatial dimensionality and temporal
parameterization of geophysical inversion is a central challenge. An
additional difficulty lies in identifying how geophysical information
should be coupled with GW models to robustly characterize subsurface
heterogeneity and hydrological processes}. In particular, it remains
unclear whether GW models can be effectively constrained by static
geophysical interpretations, or whether geophysical forward modeling
should instead be dynamically constrained by transient GW simulations.
The choice depends on several factors: (i) geological complexity,
hydrogeological context and human uses may require 3D models; (ii) the
spatial coverage and temporal resolution of the acquired datasets may
be insufficient to support fully 3D inversion; (iii) high-dimensional
models can resolve finer-scale features but require sufficient data
\citep{bai_groundwater_2021, casillas-trasvina_characterizing_2022,
neven_novel_2023}.

2D approaches provide a compromise in terms of computational efficiency
and are particularly suited to layered systems or settings where
lateral heterogeneity is limited. They are effective for understanding
processes along transects or profiles, such as SW--GW interactions or
hillslope hydrology, once the dominant flow directions and
heterogeneity patterns as well as the primary direction of change are
known or have been thoroughly determined by sensitivity analyses. In
contrast, 3D inversions provide the most comprehensive representation
of coupled transient physical processes but require substantial
computational resources and dense data coverage
\citep{linde_joint_2016}. As a result, advances in parallel computing,
workflow automation, and model reduction (see Section
\ref{sec:parameters_reduction}) or surrogate modeling approaches
\citep{asher_review_2015} are essential for applying 3D approaches to
field-scale systems with heterogeneous geology
\citep{linde_geological_2015} or with complex processes, such as
seawater intrusion \citep{steklova_joint_2017}.

An alternative approach, proposed by \citet{neven_novel_2023}, involves
iteratively coupling geological, geophysical, and groundwater FWD
models within a stochastic framework. FWD responses are computed in
lower-dimensional spaces relevant to each physical problem (the
geophysical FWD is simulated in 1D for SRT and the GW FWD in 2D) to
speed up the inversion process. This dimensional reduction exploits the
fact that different physical processes are sensitive to different
spatial dimension of the subsurface parameter field, enabling
significant computational savings with limited information loss. In
this sense, \citet{cunha_teixeira_neural_2025} for instance trains a
language-model-style neural network on many 1D FWD simulations that map
seismic velocity profiles to petrophysical and hydrogeophysical
properties. Once trained, the model can generalize these learned 1D
relationships in 3D by translating new seismic data into spatially
consistent, text-described subsurface properties, effectively lifting
1D inversion knowledge into a 3D interpretive framework. But {such
reduced}-dimensional FWD simulations cannot fully represent all 3D
processes, and the gains in computational efficiency often outweigh the
minor losses in information, particularly for large-scale or
data-limited studies. Further research into integrating FWDs of
differing dimensionalities is therefore needed to improve computational
efficiency, accuracy, and scalability in hydrogeophysical modeling,
while identifying simplifying assumptions that preserve predictive
capability \citep{schilling_beyond_2019, steklova_joint_2017}.

\subsection{Real physical integration in model vs.\ representativeness
based on the real data}

{Subsurface modeling requires a careful trade-off between physical
complexity and the representativeness of the available data.}
Increasing model complexity, such as incorporating
thermo-hydro-mechanical (THM) coupling, can improve physical realism
\citep{pandey_geothermal_2018, carcione_physics_2019,
yang_simulation_2024,haibing_characteristics_2014} 
but often limits applicability at the hillslope scale due to data
scarcity and computational demands. The coupling mechanisms of THM are
mathematically represented by parameters such as Biot coefficients and
thermal expansion or conductivity coefficients
\citep{pham_thermo-hydro-mechanical_2023,
ai_thermo-hydro-mechanical_2024, liu_simulation_2023,chambers_new_2024} 
which can be easily coupled with geophysical FWD models. A successful
application to variably saturated soils undergoing freeze--thaw cycles
demonstrates how simplifications, such as 2D shallow homogeneous
domains, isotropic deformation, and simplified hydraulic boundary
conditions that neglect lateral subsurface flow, make simulations
feasible while retaining key THM interactions
\citep{huang_numerical_2023}. The key question is how to identify and
validate simplifying assumptions that preserve predictive accuracy
while enabling practical implementation. A suitable level of model
complexity is also scale-dependent. Fully coupled THM models may be
appropriate at small scales such as the borehole application and deep
geothermal systems
\citep{hermans_geophysical_2014, shariatinik_ert_2024}, but their
applicability in the CZ remains limited due to variable saturation,
pronounced heterogeneity, high-dimensional parameter spaces, scarce
petrophysical and mechanical constraints, and pervasive uncertainty in
climate-sensitive boundary conditions
\citep{schilling_beyond_2019, neven_novel_2023, hinnell_improved_2010,
romhild_joint_2024}. In these larger domains, simpler hydro-thermal or
hydro-mechanical formulations often provide sufficient predictive skill
faced to the available data. In other words, a central challenge is to
ensure that the predictive models employed are valid for representative
elementary volumes (REV) over which the averaged equations remain
applicable and predictive. The equations that accurately reflect
reality are intrinsically linked to the REV corresponding to a given
scale. The challenge is to ensure model representativeness while
maintaining predictive accuracy through validated simplifications that
allow practical implementation. Thus, the complexity of a model should
reflect the dominant processes that can realistically be constrained at
the relevant scale.

\subsection{Managing parameter complexity and model uncertainty} 
\label{sec:parameters_reduction}

Probabilistic inversion methods are extremely computationally expensive
because they require thousands of forward simulations, especially for
complex geometries like shallow groundwater systems
\citep{ferre_being_2020}. The computational burden obviously grows with
the number of parameters to invert, making high-dimensional problems
impractical without simplification. Reducing dimensionality and
improving efficiency is essential to make probabilistic approaches
feasible for real-world hydrogeophysical applications. Therefore,
parsimonious parameterizations must also be proposed. In order to
reduce the model parameter space to be sampled, any available a priori
information should be accounted for such as the information issued from
geological models (for example, Flumy realizations). A few other
approaches to reduce computational demand have also been proposed, for
instance, the use of parallel/distributed computing technologies
\citep[e.g.][]{tavakoli_parallel_2013}; the use of efficient surrogate
models to reduce the cost of the forward problem, i.e., \mbox{reduced-order}
flow modeling \citep[e.g.][]{lu_drips_2023,gadd_surrogate_2019}, and
the reduction of the number of model parameters to infer. 

Sensitivity analysis is essential for identifying the contribution of
individual or grouped input parameters to output variance, typically
measured using Sobol sensitivity indices
\citep{sobol_sensitivity_1993}. Common indices include first-order
indices, representing isolated contributions, and total-order indices,
accounting for joint effects with other parameters
\citep{sochala_polynomial_2021}. Computational efficiency can be
improved using optimized sampling strategies
\citep{campolongo_effective_2007,saltelli_variance_2010} or
distance-based sensitivity, which is an alternative requiring less
realizations \citep{fenwick_quantifying_2014,park_dgsa_2016}. For
instance, in the poroelastic context, sophisticated MCMC
probabilistic-based variance methods have been developed, such as
E-FAST tool \citep{mesgouez_use_2017,mesgouez_characterization_2018},
to estimate the impact of the parameter perturbations on recorded
seismic data for a given experiment configuration. Those kinds of
analysis allow to reduce the number of parameters and could be also
naturally adapted to other FWDs. However, those techniques involve a
huge number of forward problem runs. An alternative approach involves
the Taguchi method, which employs statistical approaches based on
orthogonal arrays to estimate the impact of different parameters
\citep{plazolles_monte-carlo_2015,plazolles_plateforme_2017,
martin_gravity_2025} and reduce the number of parameters. This will
allow not only to reduce the computational costs of each FWD and data
inversion in terms of CPU time and disk storage but also to better
constrain the data inversions thanks to those parameter reductions and
their related impact estimates. Once the probabilistic analysis has
been done for one given reduced set of parameters associated with a
given physical domain, it is possible to study and quantify the impact
of the variations of those parameters on the reduced parameter set of
another physical domain. This also allows us to estimate the interest
of performing joint inversions or not, and to define realistic
parameter space intervals or covariance matrices to constrain
{separate,} collaborative, or simultaneous joint inversions
\citep{giraud_structural_2021,ogarko_tomofastx_2023,martin_gravity_2025}. 

In addition, adaptive reduced-order modeling approaches have been
applied successfully in \mbox{hydrogeophysics,} such as projection-based
\mbox{approaches} for resistivity imaging of solute transport
\citep{oware_geophysical_2014} and Bayesian Evidential Learning
frameworks \citep{hermans_uncertainty_2018}, for instance more recently
applied to seismic data \citep{Mreyen_full_2025}. They reduce the
parameter space by focusing only on the prediction-relevant dimensions,
rather than attempting to resolve the full subsurface model. This is
achieved through statistical learning and prior falsification, which
filters out models inconsistent with observed data. These strategies
enable efficient inversion-free prediction while preserving accuracy in
dynamic GW systems.\looseness=-1

Besides, recent advancements in machine learning, such as deep neural
networks (DNNs), have demonstrated the efficient reduction of complex
spatial patterns to manageable dimensions. DNNs simplify subsurface
models with millions of cells while maintaining geometrical complexity
\citep{laloy_trainingimage_2018, laloy_inversion_2017,
lopez-alvis_deep_2021}. 
This dimensionality reduction enables the application of gradient-based
deterministic inversion
\citep{cui_characterizing_2024, lopezalvis_geophysical_2022} or MCMC
sampling \citep{laloy_trainingimage_2018} at a reasonable computational
cost. Realistic solutions can be further constrained through
falsification approaches or the construction of assembled priors. 

Ultimately, the uncertainties can be propagated from one inverse
problem to another inverse problem \citep{gesret_propagation_2015} when
independent data are available. For example, the uncertainties on
saturation profiles could be estimated by inverting seismic data in a
first step. Then the uncertainties on hydrogeological parameters could
be estimated by inverting temperature data accounting for both the
uncertainties on temperature measurements and the uncertainties on
saturation profiles. Coupled hydrogeophysical models integrate multiple
physical processes, such as geophysical inversion and hydrogeological
simulation, which introduces multi-source uncertainty. These
uncertainties arise from measurement errors, inversion approximations,
and petrophysical relationships. Propagating errors from geophysical
data (e.g., seismic, electrical resistivity) into hydrogeological
predictions (e.g., groundwater flow, temperature distribution) is
particularly challenging because these relationships are often
nonlinear and poorly constrained \citep{linde_joint_2016,
irving_stochastic_2010}. Several strategies have been proposed to
quantify and propagate uncertainty in such coupled models. Bayesian
posterior analysis treats model outputs as random variables conditioned
on uncertain inputs. Bayesian inversion frameworks, including recent
advances using deep generative models, enable a probabilistic
characterization of uncertainty and provide credible intervals for
predictions \citep{irving_stochastic_2010, laloy_trainingimage_2018}.
Similarly, covariance propagation approaches rely on linearized error
propagation using Jacobians, which offers an efficient way to estimate
uncertainty for small perturbations, while nonlinear extensions have
been developed to better capture uncertainty amplification in complex
systems \citep{zhang_imprecise_2020, reichert_potential_2021}. 

Despite these advances, a major challenge remains: how to describe and
combine errors across parameters that differ in units, spatial
resolution, and temporal scales. Geophysical parameters such as
resistivity or seismic velocity hydrological parameters such as
hydraulic conductivity and hydrological variables such as water
pressure or saturation are not directly comparable, making unified
error metrics essential. To address this heterogeneity, several
solutions have been proposed. One approach is to normalize errors using
dimensionless metrics, such as the coefficient of variation (CV),
defined as $\mathrm{CV} = {\sigma}/{\mu}$, where $\sigma$ is the
standard deviation and $\mu$ the mean. This normalization enables
comparison across parameters with different units. Another strategy is
to represent uncertainties as probability distributions rather than raw
errors. For example, geophysical inversion may yield a posterior
distribution of resistivity, while hydrological models produce a
posterior distribution of hydraulic conductivity or Van Genuchten
parameters. These distributions can then be combined in a Bayesian
hierarchical framework, where each parameter retains its own
probabilistic space but is linked through joint likelihoods
\citep{irving_stochastic_2010, laloy_trainingimage_2018}. An
alternative approach involves transforming parameters to a common space
using petrophysical relationships, such as Archie's law, and
propagating uncertainty through these transformations using Monte Carlo
simulations or covariance-based methods, expressed as
$\mathbf{C}_{\mathrm{output}} = \mathbf{J}\,\mathbf{C}_{\mathrm{input}}
\,\mathbf{J}^T$, where $\mathbf{J}$ is the Jacobian of the
transformation \citep{tarantola_inverse_2005}. Finally, multi-objective
error metrics have been \mbox{introduced} to combine normalized errors from
both domains in a single misfit function, ensuring balanced
contributions from geophysical and hydrological data. For instance, if
geophysical inversion gives resistivity uncertainty of $\pm$10\% and
the hydrological model a hydraulic conductivity uncertainty of
$\pm$30\%, both can be normalized or represented as probability
distributions and propagated through the coupled model using Bayesian
updating or ensemble simulation.

Nevertheless, computational and practical challenges persist.
High-dimensional parameter spaces make sampling-based approaches, such
as MCMC or ensemble methods, computationally expensive. Nonlinear
petrophysical relationships amplify uncertainty during propagation, and
visualization of uncertainty remains difficult, requiring improved
tools for probability mapping, uncertainty envelopes, and multi-scale
representation \citep{hermans_uncertainty_2018}. Future research should
focus on developing unified frameworks that can handle multi-scale,
multi-physics uncertainty propagation
\citep{gomesgoncalves_uncertainty_2025}.

\subsection{Model-guided monitoring design}

Effective integration of geophysical and hydrological data into
groundwater models requires not only robust coupling approaches but
also strategic monitoring design to optimize observation types,
locations, and frequencies. However, identifying the most informative
dataset under uncertainty remains a fundamental challenge
\citep{gates_worth_1974,freeze_hydrogeological_1992}. 

Data-Worth Analysis (DWA) provides a framework to quantify the value of
geophysical data in improving hydrological model predictions. By
comparing model outcomes with and without specific geophysical
datasets, DWA helps assess their influence on parameter estimation and
forecast reliability. Value-of-Information frameworks extend DWA by
explicitly incorporating decision-theoretic and economic considerations
into monitoring design
\citep{alfonso_coupling_2012,trainor-guitton_geophysical_2014}.
\citet{feyen_framework_2005} framed data worth within Bayesian decision
analysis for groundwater management, illustrating how economic
objectives determine whether additional hydraulic conductivity
measurements justify their collection costs.

Multi-objective monitoring network design approaches combine data
assimilation with optimization algorithms to identify observation
locations and types that maximize information gain while respecting
budget constraints. \citet{thibaut_comparing_2022} applied Bayesian
Evidential Learning (BEL) within an optimal experimental design
framework to compare the value of borehole temperature logs versus
time-lapse electrical resistivity tomography (ERT) for monitoring 4D
thermal plumes, addressing the computational challenges and tradeoffs
between invasive and non-invasive monitoring methods while identifying
optimal sensor combinations that minimize prediction uncertainty under
budget constraints. \citet{durney_multipurpose_2025} applied
multi-purpose DWA to a coupled SWAT-MODFLOW-RT3D model, evaluating the
spatial worth of nitrate concentration measurements and SkyTEM-derived
geophysical data for improving predictions of both water quantity and
quality. Similarly, \citet{vilhelmsen_extending_2018} extended the
approach to multi-objective contexts, allowing the selection of
observations that support several forecasting goals simultaneously.

Adaptive sampling strategies, such as those used in time-lapse ERT
surveys \citep{wilkinson_adaptive_2015}, adjust measurement
configurations dynamically based on prior observations. Surrogate-based
optimization approaches, including proper orthogonal decomposition and
Gaussian process modeling, allow for efficient evaluation of monitoring
scenarios in computationally intensive models
\citep{gosses_robust_2021}. \citet{mern_intelligent_2023} proposed
sequential decision-making approaches using Monte Carlo planning
methods for mineral exploration, demonstrating that closed-loop
sequential information gathering is superior to non-sequential schemes
while explicitly addressing the combinatorial explosion problem and
computational burden inherent in optimal experimental design. Most
data-worth studies focus on single prediction targets (e.g.,
contaminant arrival time), while water resources management
increasingly requires multi-objective predictions (water quantity,
quality, ecosystem services) that may demand different optimal
monitoring strategies {\citet{vilhelmsen_extending_2018}}. Future
research should focus on refining these approaches to ensure that
monitoring efforts effectively enhance model reliability while
remaining feasible under financial and logistical constraints.

\subsection{From hillslope models to catchment-scale simulations: the
role of change-of-support}

Simulating water resources, especially on a watershed scale,
necessitates the integration of datasets with various spatial supports,
including simulated water and heat fluxes at {hillslope} model grid
resolutions, coarse-scale remote sensing observations (e.g., vegetation
indices, soil moisture, topography) \citep{lubczynski_remote_2024},
fine-scale geophysical interpretations (ERT, GPR, seismic methods), and
point-scale records such as piezometric levels and temperature.

Except recent advances showing promising results at the regional scale
with ambiant noise seismology and geodesy
\citep{denolle_ambient_2025,lu_mapping_2025}, ERT and {active} seismic
methods are not usable at the full catchment scale. However, their
hillslope-scale resolution remains still beneficial when used with
scalable geophysical methods and change-of-support techniques.
Fine-scale information derived from ERT and SRT can constrain or
parameterize larger-scale Frequency-Domain Electromagnetic (FDEM) or
Time-Domain Electromagnetic (TDEM) surveys, by providing information
about lithological contrasts, hydrological connectivity, and shallow
structural patterns that EM methods alone cannot resolve. Furthermore,
groundwater models utilizing high-resolution geophysical data can
estimate stream--aquifer interactions and aquifer recharge on larger
scales, thereby enhancing the physical realism of catchment-scale
simulations.

To move from hillslope-scale simulations, such as water fluxes or
SW--GW fluxes, to watershed domains, or to integrate interpretations
from remote sensing and hydrogeophysics, change-of-support approaches
are crucial
\citep{emerychangesupport2009,zaytsevchangesupport2016}. These
methods help reconcile differences in data support, ensuring consistent
interpretation of variance, correlation, and uncertainty at the larger
scale \citep{de_marsily_dealing_2005}. Methods such as block kriging,
conditional simulation, and nonlinear geostatistical approaches (e.g.,
disjunctive kriging, uniform conditioning) help bridge between
hillslope and watershed scales
\citep{demougeot-renard_geostatistical_2004}. Conditional simulation
involves generating multiple realizations of the random field
representing the target attribute, followed by regularization to block
support \citep{zaytsevchangesupport2016}. While powerful, this
method is computationally demanding, requiring specification of
finite-dimensional distributions and conditioning to sparse data.
\citet{benoit_stochastic_2021} addressed the change-of-support problem
by up-scaling local quasi-point conductivity fields to block-scale
conductivity tensors for GW flow modeling. The use of nonlinear
geostatistical methods, such as disjunctive kriging or uniform
conditioning provide an analytical solution to the change-of-support
problem \citep{emerychangesupport2009}.
\citetalias{emerychangesupport2009} and
\citet{zaytsevchangesupport2016} describe two generalizations of
{Discrete Gaussian Models} (DGMs) for irregular grids. DGMs ensure
accurate reproduction of block-to-block covariance and marginal
distributions and simplify computation. These models are particularly
relevant for unstructured grids, such as those used in reservoir
modeling, and are compatible with central limit theorem behavior under
variance reduction. Additional research is needed to explore a wider
range of physical processes and to derive scaling laws for the
occurrence and characteristics of GW flow systems at various scales.

\section{Conclusions}

Advancing hydrogeophysics requires bridging the gap between its
theoretical potential and its operational application to support
sustainable critical zone science, GW management, resilience under
global change and geothermal applications. One of the primary
objectives of hydrogeophysics is to develop monitoring and forecasting
tools for resolving the spatiotemporal variability of shallow GW under
global change. This requires resolving the transient physical processes
that govern water and heat fluxes, vectors of solute transport and
transformation within the critical zone. This paper details three
different geophysical methods commonly used at the hillscope scale.
Many tasks remain to be done to improve the imaging capabilities in
heterogeneous transient contexts regarding noise analysis and data
integration. This paper underscores a paradigm shift from static to
transient inversion, where coupled and joint approaches provide new
opportunities to resolve subsurface heterogeneity and GW variability.

Progress will depend on addressing parameter interdependence,
petrophysical uncertainty, and mismatches in the spatial and temporal
scales of \mbox{multi-physics} datasets. Geophysical data play a critical role
in improving GW model initialization and boundary conditions by
providing spatially distributed information that is otherwise difficult
to obtain through direct measurements. A central challenge lies in
determining the appropriate dimensionality for geophysical inversion
and its coupling with GW models. In some cases, 1D or 2D geophysical
inversions can be better constrained by using state variables simulated
in 3D GW models, thereby linking scales and enhancing the
representation of transient processes. However, such workflows must be
carefully validated to ensure that the geophysical model remains
consistent with hydrological reality and that unnecessary 3D inversion
can be avoided. 

Data availability remains limited by budgetary, logistical, and
technical constraints, which often necessitate conceptual
simplifications such as reducing the number of parameters,
dimensionality, or the detail in initial and boundary conditions.
Inversion dimensionality should therefore be adapted to site-specific
contexts, balancing complexity with representativeness, while
maintaining flexible rather than prescriptive workflows. Such
simplifications should not be adopted as default assumptions but rather
as deliberate, evidence-based choices supported by experiments,
sensitivity analyses, or mechanistic modeling. The principle of
parsimony remains valid, provided it does not compromise predictive
capability. 

At the same time, models should guide the design of experimental
networks and measurement strategies, ensuring that data collection
targets the most sensitive and uncertain processes. Geostatistical
frameworks such as SPDE and change-of-support approaches require
further research to strengthen hydrogeophysical inversion workflows,
not only for integrating multi-source data and interpolating poorly
constrained initial and boundary conditions, but also for using outputs
from 3D GW models to better constrain and inform 2D geophysical
inversions or to interpolate between multiple 2D inversions within a 3D
hydrological framework. Addressing stakeholder needs will require
extending model to the watershed scale, while integrating insights
gained at the hillslope scale and leveraging remote sensing. Future
research on change-of-support remains essential.

Meeting these challenges will then depend on scalable inversion
strategies that combine inversion, probabilistic frameworks, and
machine learning surrogates. Field observatories linked to open-access
databases and equipped with advanced monitoring technologies such as
CZNet, OZCAR, and TERENO \citep{gaillardet_ozcar_2018} are also
essential to (i) build mechanistic models grounded in observations,
(ii) test and validate sensors and methodologies, and (iii) identify
which processes truly require high-resolution space--time data. By
reducing uncertainty, optimizing dimensionality, and embedding field
validation, hydrogeophysics can become a cornerstone of critical zone
science and a key enabler of resilience in the Anthropocene.

\section*{List of abbreviations}

\begin{inftab}
\fontsize{10}{11.8}\selectfont
\begin{tabular}{ll}
\tbody
CZ & Critical xone \\
GW & Groundwater \\
SW & Surface water \\
SW--GW & Surface water--groundwater interactions \\
VZ & Vadose zone \\
BC & Boundary conditions \\
ERT & Electrical resistivity tomography \\
SRT & Seismic refraction tomography \\
WRC & Water retention curve \\
FWD & Forward model \\
$V_P$ & P-wave velocity \\
$V_S$ & S-wave velocity \\
MASW & Multichannel analysis of surface waves \\
MCMC & Markov chain Monte Carlo\\
SHI & Sequential hydrogeophysical inversion \\
CHI & Coupled hydrogeophysical inversion \\ 
JHI & Joint hydrogeophysical inversion \\
PBI & Process-based inversion \\
DNNs & Deep neural networks \\
MPS & Multiple-point simulation \\
SPDE & Stochastic partial differential equations \\
PGS & Plurigaussian geostatistical simulation \\
DGM & Discrete Gaussian model\\
BSS & Bayesian sequential simulation method \\
DWA & Data-worth analysis \\
1D & One-dimensional (space) \\
2D & Two-dimensional (space) \\
3D & Three-dimensional (space) \\
4D & Four-dimensional (space and time)
\end{tabular}
\end{inftab}

\section*{Acknowledgments}

The manuscript was written through contributions of all authors.
The authors acknowledge the support of the French Agence Nationale de
la Recherche (ANR) under grant ANR-23-CE01-0004 (project GWSBOUND).
They also acknowledge the infrastructure and resources provided by
OZCAR-RI (Observatoires de la Zone Critique -- Recherche et
Infrastructure). The authors acknowledge the support of Institutes UMR
7619 METIS, UAR3455 ECCE TERRA, Sorbonne Universit\'e, CNRS and Mines
Paris - PSL (PSL University).

\CDRGrant[ANR]{ANR-23-CE01-0004}

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations.

\back{}

\appendix{}

\def\appendixlabel{Appendix \Alph{section}}

\section{Geological constraints on GW fluxes: example of sediment
heterogeneity in alluvial aquifers}\label{app:geology}
\subsection{Structure and heterogeneity of alluvial aquifers in
north-western Europe}

River dynamics lead to the sediment deposition and formation of
sedimentary bodies of varying lithologies and geometries, which
determine the sedimentary heterogeneity of alluvial plain environments.
These sedimentary heterogeneities, which correspond to contrasts in
lithological and grain-size properties
\citep{koltermann_heterogeneity_1996}, also result in contrasts in
hydraulic parameters (e.g.\ hydraulic conductivity, porosity) and
transport parameters (e.g.\ effective diffusion coefficient, retardation
factor, diffusivity, thermal conductivity, heat capacity). This has led
to the concept of hydrofacies, developed within the hydrogeology
community \citep{poeter_influence_1990}. Sediment heterogeneities
therefore influence preferential flow paths and residence times of
water in the alluvial plain, as well as on GW-river exchanges
\citep{fogg_groundwater_1986, anderson_sedimentology_1999,
webb_simulation_1996, hornung_reservoir_1999, klingbeil_relating_1999,
bianchi_spatial_2011, heinz_heterogeneity_2003,
fleckenstein_river-aquifer_2006,miall_geology_1996,stonedahl_multiscale_2010,flipo_continental_2014}. 
The streambed heterogeneities, coupled with the longitudinal variation
of the bed, impact the dynamics of the SW--GW exchanges by creating
flow recirculation \citep{cardenas_impact_2004}. Sediment
heterogeneities are the result of variability in time and space in the
functioning of the fluvial systems, and more specifically in the energy
of the depositional environment. They exist and can be described on
different space and time scales, from lamina resulting from turbulent
processes on the scale of seconds, to the basin controlled by climate
and tectonics, which modify flow and sediment supply conditions on a
regional scale \citep{miall_geology_1996}.\looseness=1

Most of the alluvial plains of north-western Europe were shaped by the
alternating glacial and interglacial periods that characterized the
Late Quaternary from the Middle Pleistocene
\citep{antoine_chronostratigraphy_2010}. These alternations are
reflected in the response of fluvial systems to incision or
aggradation. It is particularly during periods of climatic transition
that fluvial systems become unstable and undergo profound morphological
transformations. Climatic optima and glacial maxima, on the other hand,
are periods of sediment reworking. In coastal settings, cutting and
filling respond strongly to base-level fluctuations driven by
glacioeustatic sea-level changes \citep{schumm_river_1993,
dalrymple_incised_2006, cardinal_interplay_2024}. Upstream of the limit
of eustatic influence, changes in vegetation cover and hydrological
functioning essentially modulate the ratio of water to solid discharge,
and thus the timing of incision or aggradation of the fluvial system.
Generally, large incisions are observed at interglacial-glacial
transitions, followed by strong aggradation during the glacial period
by high-energy braided fluvial systems. A more moderate incision marks
the transition to the interglacial, where lower-energy systems, such as
meandering channels, rework the glacial braided deposits
\citep{vandenberghe_relation_2002}.

While braided river systems can temporarily store large volumes of fine
sediments during low-water periods, these are rapidly remobilized
during floods \citep{navratil_assessment_2010}. As a result, their
preservation potential is relatively limited. The sediment
heterogeneity of alluvial deposits formed by these systems results
essentially from the degree of mixing between sands and gravels that
occurs in unit and compound bars, from open-framework gravel beds to
fine sand drapes \citep{bluck_structure_1979, bridge_fluvial_2006,
klingbeil_relating_1999, heinz_heterogeneity_2003}. 

Conversely, meandering fluvial systems generate highly heterogeneous
deposits, with contrasting lithological properties expressed at
different spatial scales \citep{gouw_alluvial_2007}. At the kilometer
scale, fine overbank deposits and abandoned channels filled with high
organic content clays produce sharp lithological contrasts with the
coarser channelized facies (Figure~\ref{fig:hetero_sedi} Level 1 and
2). The low width-to-depth ratio of meandering channels, compared to
braided channels, allows the deep reworking of braided alluvial
deposits associated with older glacial periods. At the scale of several
tens to hundreds of meters (Figure~\ref{fig:hetero_sedi} Level 3), the
internal structures of sand and gravel bars reflect the distribution of
sediment load in the water column, the succession of different flow
stages, and the morphodynamic interactions with other bars and
cross-bar channels \citep{bridge_fluvial_2006}. Their sizes, highly
variable but proportional to the channel size, range from several tens
to several hundreds of meters in width, and in height reflect the
maximum channel depth. Grain segregation in bedload and turbulent
fluctuations of the flow produce heterogeneous cross-stratification
within bedforms at the centimeter scale
\citep{allen_classification_1963, allen_bed_1966}
(Figure~\ref{fig:hetero_sedi} Level~4). 

\begin{figure*}
\includegraphics{fig07}
\caption{\label{fig:hetero_sedi}Hierarchical levels of sedimentary
heterogeneity at aquifer-river interfaces (modified from 
\citet{jordan_hierarchical_1992}).} 
\end{figure*}

\subsection{Field characterisation and numerical modeling of alluvial
heterogeneity}

The heterogeneities of alluvial aquifers can be assessed directly in
the field through outcrops, trenches, or boreholes
\citep{heinz_heterogeneity_2003, hornung_reservoir_1999,
klingbeil_relating_1999, bertran_late_2025, malenda_grain_2019},
or indirectly via geophysical imaging methods. While sedimentological
analysis provides detailed lithofacies descriptions along 1D or 2D
profiles, the data are localized and challenging to interpolate.
Geophysical methods, on the other hand, deliver broader spatial 2D/3D
data but often lack resolution for smaller-scale heterogeneities.

Numerical methods offer the possibility to create fully spatialized
models of alluvial aquifers \citep{bayer_three-dimensional_2011,
comunian_three-dimensional_2011},
resolving the different scales of sediment heterogeneities, and
potentially conditioned by field data
\citep{zappa_modeling_2006, bianchi_lithofacies_2016,
sun_characterization_2008, teles_comparison_2004, zhang_impact_2013}. 
Geostatistical methods reproduce the geometric characteristics of an
alluvial reservoir, without taking into account the reservoir
construction history or the processes behind it. In object-based
methods, each sedimentary body is treated as an \mbox{object,} discretized on
a grid or continuously modeled by surfaces
\citep{jacquemyn_surface-based_2019}. The domain is filled using
probability laws describing both the spatial distribution and the size
and orientation of the objects, together with laws of
attraction/repulsion between objects. In pixel-based methods, training
images are used to define the geostatistical characteristics of
sedimentary facies. Realizations are then performed following facies
ordering rules \citep{beucher_truncated_2016,tahmasebi_multiple_2018}.
Geostatistical models can be easily conditioned by field data, and low
computation times enable a large number of realizations to be
generated. On the other hand, results sometimes lack geological
realism. Process-based or genetic models aim to reproduce the evolution
of a fluvial systems and the associated sediment deposits through a
simplification of the physical laws of flow and sediment transport
\citep{camporeale_hierarchy_2007}. Because computation times can be
very long, deterministic models are generally limited to simulating the
morphodynamic evolution of the active channel of a fluvial system on
time scales of a few hundred to thousands of years
\citep{van_de_lageweg_preservation_2016}. Finally, mixed
stochastic-genetic models (forward
model---\xcitealp{lopez_modelisation_2003}{2003} or
reverse---\xcitealp{parquer_combined_2020}{2020}) use a deterministic approach
to compute the channel-form evolution through time, but do only mimic
the sedimentary processes with geometrical rules to build the alluvial
plain heterogeneity.

\subsection{Impact of alluvial plain heterogeneity on ground water
fluxes---a synthetic case example}

Meandering fluvial systems are characterized by highly sinuous channels
(sinuosity ${>}$ 1.5). The evolution of meander loops in time can lead to
neck cutoffs \citep{gao_exploring_2024} and abandoned oxbow lakes that
are filled with fine sediments and organic matter
\citep{jordan_hierarchical_1992}. Head gradients along meander loops
induce significant SW--GW fluxes across pointbars
\citep{cardenas_effect_2008, cardenas_streamaquifer_2009,
boano_sinuosity-driven_2006, boano_linear_2010}. 
On the contrary, mud plugs in abandoned meander loops create extensive
permeability barriers within the alluvial plain. Here, we use the Flumy
genetic sedimentary model \citet{lopez_processbased_2008},
\url{https://flumy.minesparis.psl.eu/}) \mbox{coupled} with hydrogeological
model to investigate with synthetic cases the effects of alluvial plain
heterogeneity on GW flow and water level in an alluvial aquifer.

Flumy is a model that uses stochastic and process-based approaches to
simulate meandering channelized reservoirs. It describes channel
evolution through time using physical equations from hydraulic studies
\citep{ikeda_bend_1981}, with constant section, slope, and flow.
Sedimentary facies are deposited at specific locations in the
floodplain, depending on the distance to the channel or a specific
process. Sandy point-bar is deposited inside meander loops, while
floods produce fine-grained overbank deposits, contributing to channel
aggradation. Levee breaches are at the origin of crevasse splays, and
channels abandoned after neck cutoff or avulsion are filled with
mudplug facies.

Here, Flumy generates a simple 7 km-long, 3 km-wide, 10 m-thick
floodplain heterogeneity model. The initial domain is filled with sand
facies. The 50 m-wide and 6 m-deep channel is initially straight with
small curvature perturbations and flows from west to east. Flumy is run
without overbank flood (no vertical accretion) nor avulsion for a
certain number of iterations, until the sinuosity is well developed
($D_{\mathrm{curve}}/D_{\mathrm{straight}}=2.2$), producing meander
neck cut-offs and mud plugs. The lateral migration of the channel
meanders leads to the deposition of sandy point bars, a facies
identical to the initial domain fill. Thus, mud plugs in meanders
abandoned after neck cut-off are the only source of heterogeneity in
the floodplain. From this heterogeneous model
(Figure~\ref{fig:sed_setup}d), a reference homogeneous model is
\mbox{derived} by replacing the mud plug facies by sand
(Figure~\ref{fig:sed_setup}a).\looseness=1

The hydrofacies values are derived from the literature
\citep{morris_summary_1967, urumovic_effective_2014}. The pointbar sand
has a hydraulic conductivity of $2.9\times 10^{-3}$~m${\cdot}$s$^{-1}$
while the mudplug has a conductivity of $1.6\times
10^{-9}$~m${\cdot}$s$^{-1}$.

The GW flow simulations for both models, homogeneous and heterogeneous,
are performed using the 8.1 version of Feflow
\citep{diersch_feflow_2013}. The aquifer is supposed to be in a
confined state. The simulation domain is subjected to a 
$0.4\permil{}$ regional hydraulic gradient with constant
Dirichlet conditions of 9.5 and 7.5~m upstream and downstream. The
river's hydraulic head decreases linearly along the river channel from
9~m to  7~m, representing that it is draining the aquifer. A 30-day
flood wave scenario was applied to the riverline BC, with a 3~m
increase in the first 15~days and a return to initial state level in
the last 15~days.

The resulting simulated head and Darcy flux fields in
Figure~\ref{fig:sed_setup} for the peak flow condition ($t=15$ days)
are shown for a location between two active meanders 
(Figure~\ref{fig:sed_setup}d), with the occurrence of a mudplug filled
meander.

These results, below the mudplug (Figure~\ref{fig:sed_setup}b,~c),
and cross-section view (Figure~\ref{fig:sed_setup}e,~f), exhibit a
doubling of Darcy velocity (from 10 to 20~cm/day) below the mudplug,
and a modification of the flow direction. The heterogeneity affects the
direction of GW flow beneath the mudplug: the mudplug's barrier effect
tends to reorient the flow in the direction that is orthogonal to the
mudplug, pushing downward flow below it. These results, in terms of
flow direction and magnitude, confirm that the occurrence of large
sedimentary heterogeneity, such as abandoned meander mud plug, can
generate a preferential flow path for GW flow. The influence of
significant heterogeneity on hydraulic head is restricted to the
heterogeneity perimeter, however its effect on flow velocity extend
across the full inter-meander area. This means that it could be ``easily
missed'' by a classical data set where only a few piezometers would be
available, highlighting the need for geophysical data. This implies
that such features could be overlooked in conventional data with sparse
piezometer coverage, underscoring the importance of incorporating
geophysical data.

The common approach for GW modeling involves calibrating hydrodynamical
parameters by comparing simulated with observed hydraulic head,
neglecting flow velocity. Consequently, the influence of sedimentary
heterogeneity on flow velocity may pose a significant challenge for
applications like contamination and heat transport or flow
quantification.

\section{Physical equations}\label{app:equations}
\subsection{GW flow equations}\label{app:flow_eq}

The mass balance equation used to calculate GW (GW) flow is: 
{\begin{equation}\label{eq:flow_eq}
\nabla \cdot \left(\dfrac{k \bar{k}_r(S_w) \rho _w}{\mu _w} (\nabla p -
\rho _w g) \right) = \dfrac{\partial (\rho _w n S_{w})}{\partial t} + 
\rho _w q,
\end{equation}}\unskip
where: $t$ is time (s), $p$ is the pressure (Pa), $S_{w}$ is the water
saturation,  $k$ is the intrinsic permeability (m\tsup{2}), 
$k_{r}(S_w)$ is the relative unsaturated hydraulic conductivity, given
by the water retention curve, $g$ is the acceleration due to gravity
(m${\cdot}$s$^{-2}$), $\mu_w$ is the dynamic viscosity of water
(Pa${\cdot}$s),  $\rho_{w}$ is the density of water (kg${\cdot}$m$^{-3}$),
$n(x, y, z)$ is porosity (-), $q$ is a volumetric water source or sink
(s$^{-1}$). 

\subsection{Heat transport equations} \label{app:heat_eq}

The heat transport is described by the advection-diffusion equation: 
{\begin{equation}\label{eq:heat_transport} 
\nabla [\lambda\,\nabla T -\rho_{w}\,C_{w}\,q\,T] = \rho\, C\,
\dfrac{\partial T}{\partial t},
\end{equation}}\unskip
where $T$ is temperature ($K$), $\lambda$ is the thermal conductivity
of the porous medium (water, air $+$\ solid)
(W${\cdot}$m$^{-1}{\cdot}$K$^{-1}$), $\rho$ is the density of the porous
medium (kg${\cdot}$m$^{-3}$), $C_{w}$ is the specific heat capacity of
water (J${\cdot}$kg$^{-1}{\cdot}$K$^{-1}$), $C$ is the specific heat
capacity of the porous medium (J${\cdot}$kg$^{-1}{\cdot}$K$^{-1}$) and
$q$ is the specific discharge (m${\cdot}$s$^{-1}$). 

\subsection{Electrical equations}\label{app:ert}

The equations and theoretical models describing electrical current flow
in materials are detailed in \citet{telford_applied_1990}, particularly
for applications in electrical resistivity sounding and tomography.
These models combine Ohm's law with the governing equation for
electrical potential, resulting in Laplace's or Poisson's equations.

Ohm's law in a conductive medium is expressed as:
{\begin{equation}\label{eq:ohm}
\mathbf{J} = \sigma \mathbf{E},
\end{equation}}\unskip
where $\mathbf{J}$ is the current density (A/m$^2$), $\sigma$ the
electrical conductivity (S/m) wich is the inverse of the electrical
resistivity, and $\mathbf{E}$ the electric field (V/m).

The electric field relates to the electrical potential $\phi$ by:
{\begin{equation}\label{eq:efield}
\mathbf{E} = - \nabla \phi.
\end{equation}}\unskip

Combining these and applying the principle of conservation of charge
(no accumulation of charge in steady state):
{\begin{equation}\label{eq:charge}
\nabla \cdot \mathbf{J} = Q,
\end{equation}}\unskip
where $Q$ is the source term (current density per volume).

Substituting Equations (\ref{eq:ohm}) and (\ref{eq:efield}) into
(\ref{eq:charge}), we get the governing equation for the electrical
potential:
{\begin{equation}\label{eq:governing}
\nabla \cdot (\sigma \nabla \phi) = -Q.
\end{equation}}\unskip

In regions without current sources, $Q = 0$, this reduces to Laplace's
equation:
{\begin{equation}\label{eq:laplace}
\nabla \cdot (\sigma \nabla \phi) = 0.
\end{equation}}\unskip

\subsection{Poro-elastic wave equations}\label{app:poroelastic_eq}

Porous materials are made of a solid phase (called the frame) and of a
fluid phase, and can be considered as an interconnected network of
pores inside the solid \citep{pride_seismic_2004}. 

Air and/or water can play the role of the fluid and grains the solid.
Poroelastic materials are most of the time modeled using the Biot
theory \citep{biot_theory_1956,biot_theory_1956-2}. The differential
or ``strong'' formulation of the poroelastic wave equations can be
written as in \citet{carcione_wave_2007} and \citet{carcione_wave_2014}. 

The related system of equations has seven wave eigenvalues related to
seven wave velocity modes (instead of five for the elastic case). Those
wave velocities are $\pm V_{p\mathrm{FAST}}$, $\pm V_{p\mathrm{SLOW}}$,
$\pm V_{s}$ and 0. The shear velocity $V_s$ and the fast and slow
P-wave velocities ($V_{p\mathrm{FAST}}$ and $V_{p\mathrm{SLOW}}$) can
be expressed as \citep{sidler_pseudospectral_2013}:
{\begin{eqnarray}
V_s &=& \sqrt{\dfrac{\mu}{a_1}}, \Seqnsplit
V_{p_{\mathrm{FAST}}} &=& \sqrt{\dfrac{-b_1 + \sqrt{\Delta}}{2 a_1}}, \Seqnsplit
V_{p_{\mathrm{SLOW}}} &=& \sqrt{\dfrac{-b_1 - \sqrt{\Delta}}{2 a_1}},
\end{eqnarray}}\unskip
where 
{\begin{eqnarray*}
a_1 &=& \rho_{11} \rho_{22} - \rho_{12}^2, \\
b_1 &=& -S \rho_{22} - R \rho_{11} + 2 ga \rho_{11}, \\ 
c_1 &=& S R - ga^2, \\
\rho_{11} &=& \rho + \rho_f n (a - 2), \\
\rho_{12} &=& n \rho_f (1-a), \\
\rho_{22} &=& a \rho_f n, \\
S &=& \lambda + 2 \mu, \\
R &=& M n^2, \\
\Delta &=& b_1^2 - 4 a_1 c_1, \\
ga &=& M n (\alpha - n). 
\end{eqnarray*}}\unskip

The different parameters of the poroelastic model are: $n$ (porosity),
$\rho$ (density of the satured medium), $\rho_s$ (density of the solid),
$\rho_f$ (density of the fluid), $\rho_w$ (apparent density), $K_{s}$
(incompressibility modulus of the solid matrix), $K_{f}$
(incompressibility modulus of the fluid), $K_{fr}$ (incompressibility
modulus of the porous frame), $\kappa$ (permeability of the solid
matrix), $\eta$ (fluid viscosity), $a$ (tortuosity), $\lambda$ (Lam\'e
coefficient of the saturated matrix), $\lambda_s$ (Lam\'e coefficient
in the solid matrix), $\mu$ (shear modulus). Some of these parameters
are a combination of the others.

Biot's characteristic frequency $f_c$ defines the transition between
two poroelastic regimes (with or without attenuation) and is given as
follows (see \citet{biot_theory_1956-2}, \citet{carcione_wave_2007} and
\citet{morency_spectral-element_2008}):
{\begin{equation}\label{fc_regime}
f_c = \min \left(\dfrac{\eta n}{2 \uppi a \rho _f \kappa}\right).
\end{equation}}\unskip

The maximum frequency range $f_{\max}$ of the source is such that
$f_{\max}<f_c$. Therefore, in the experimental and numerical modeling
of unconsolidated granular media under study, we choose to stay in the
poroelastic regime without attenuation. 

\printbibliography
\refinput{crgeos20250034-reference.tex}

\end{document}
