Combining Models of Root-Zone Hydrology and Geoelectrical Measurements: Recent Advances and Future Prospects

Recent advances in measuring and modeling root water uptake along with refined electrical petrophysical models may help fill the existing gap in hydrological root model parametrization. In this paper, we discuss the choices to be made to combine root-zone hydrology and geoelectrical data with the aim of characterizing the active root zone. For each model and observation type we discuss sources of uncertainty and how they are commonly addressed in a stochastic inversion framework. We point out different degrees of integration in the existing hydrogeophysical approaches to parametrize models of root-zone hydrology. This paper aims at giving emphasis to stochastic approaches, in particular to Data Assimilation (DA) schemes, that are generally identified as the best way to combine geoelectrical data with Root Water Uptake (RWU) models. In addition, the study points out a more suitable objective function taken from the optimal transport theory that better captures complex geometry of root systems. Another pathway for improvement of geoelectrical data integration into RWU models using DA relies on the use of stem based methods as a leverage to introduce more extensive root knowledge into RWU macroscopic hydrological models.


FUNDAMENTAL CONCEPTS, ISSUES, AND OPEN QUESTIONS
The benefits of combining hydrological modeling with abundant information derived from different observations have been successfully demonstrated at least over the past 20 years. In order to support hydrological models of soil-plant interactions, observations are generally retrieved from plants (e.g., water potential measurements- Sus et al., 2014), from the atmosphere (e.g., remote sensing measurements- Gelsinari et al., 2020;Lei et al., 2020), and from the soil (e.g., soil moisture- Han et al., 2012;Gruber et al., 2015;Hu et al., 2017Hu et al., , 2019Liu et al., 2017;Zhang et al., 2018). However, most of the studies neglect the root architectural traits due to the difficulty in collecting the relevant data. This is a major issue considering the fundamental role of roots in water regulation in the Earth Critical Zone (ECZ) (Newman et al., 2006;Srayeddin and Doussan, 2009). Existing studies combining observation of roots and root-system models are limited to laboratory experiments via X-ray or Magnetic Resonance Imaging (MRI; Koch et al., 2019), but this is not the scale suitable for field and real case applications.
The use of geophysical methods has been long advocated as a possible approach to image roots and root water uptake (e.g., Werban et al., 2008). In particular, geoelectrical methods offer a reliable tool to measure soil-roots interactions. These methods rely on Electrical Resistivity (ER), which is very sensitive to soil moisture (Michot et al., 2003) and salinity content (Kemna et al., 2002). In particular, time lapse measurements have demonstrated their effectiveness in characterizing processes inherent to water movement due to root water uptake (e.g., Garré et al., 2012). Alone, geoelectrical measurements are limited to the observation of state variable changes (moisture content, solute concentration, temperature). This requires that a link with the hydrological processes inducing these state variable changes be made in order to obtain useful information from the data. Also, advances in characterizing the root zone directly using stem-based methods (capacitive and Mise-à-la-masse) Peruzzo et al., 2020) are potentially a useful tool in refining hydrological models of roots (Moreno et al., 2015;Cassiani et al., 2020;Cimpoiaşu et al., 2020;Garré et al., 2021).
Great progress has been made over the past decade to combine geophysical observation data with hydrological models (Camporese et al., 2010(Camporese et al., , 2011Tso et al., 2020). Relatively few hydro-geophysical applications, though, have been focused on root system parametrization Cassiani et al., 2016;Mary et al., 2020;Rao et al., 2020). In a series of studies, (Kuhl et al., 2018), Kuhl et al. (2021) propose a novel coupled hydrogeophysical approach to estimate evapotranspiration (ET) and root depth from ER, and soil water content combining both soil water content from Time Domain Reflectometry (TDR) and ER measurements with a hydrological and plant model. The data integration using a coupled inversion showed the advantages of ER tomography for extrapolation beyond the local scale of soil water content sensors.
Despite these advancements, many uncertainties due to model errors, data errors, inversion errors, transformation from geophysical to hydrological quantities, and initial conditions need to be accounted for (Rubin and Hubbard, 2005). Stochastic methods provide a systematic framework for assessing or handling some of the complexities that arise in fusing disparate data sets. Their ability to include structural and parametric error distributions make them particularly attractive for application to the problem of dynamic parameter estimation (Linde et al., 2017;Tso et al., 2020).
The inevitable difficulties along the way do not change the value of these approaches for blending geophysical data with hydrological models, as: (i) geophysical data offer spatially extensive information, that can be also collected frequently in time; (ii) the "calibration" of hydrological models identify their own governing parameters on the basis of the (time-lapse) geophysical data, thus giving sense to the signals recorded by the geophysical measurements. Geophysical time lapse survey focused on ways to extract relevant hydrological information using time series analysis techniques (Binley et al., 2015). Relevant time scales for root hydrology range from single measurement to daily, seasonal or even continuous monitoring.
This review paper focuses on the combination of geoelectrical observations and hydrological modeling, with the aim of characterizing the active root zone. We discuss the approaches of model parametrization according to their degree of integration from simple trials to fully coupled inversion. An emphasis is given to Data Assimilation (DA) schemes as a stochastic estimation problem. Finally, the review paves possible ways of improving hydrogeophysical models based on the current state of the art of geoelectrical methods view data assimilation.

Root Water Uptake Models
Many authors underlined the importance of identifying the best approach to describe RWU. Given the resolution of geoelectrical methods, macroscopic models (e.g., Muma et al., 2013;Vanderborght et al., 2021) describing RWU are often preferred to microscopic (or functional structural) ones (e.g., Javaux et al., 2008). These two types of models are usually based on solving a modified form of the 3D Richards equation to account for a certain root distribution and root/trunk xylem/stomatal hydraulic conductance (Volpe et al., 2013). The two approaches differ in their complexity in considering the inclusion of explicit root hydraulic features to simulate water flow toward individual roots (Doussan, 1998;Javaux et al., 2008). Currently, the link between hydraulic conductance and electrical properties of individual roots is not sufficiently known, limiting the use of combining geoelectrical measurements with microscopic models. In macroscopic models, RWU intervenes as a sink term in the soil water flow equation without solving for flow toward individual roots and thus the root system can be considered using a density formulation from empirical models (Vrugt et al., 2001). Examples of successful applications of macroscopic models are described by Couvreur et al. (2012), Manoli et al. (2014) and Camporese et al. (2015b).
Typical application of macroscopic RWU models assumed the root-zone averaged saturation (Manoli et al., 2014), making the hydrological model segments' length comparable with geophysical ones. Note that one way to simplify the description of microscopic water flow is to upscale hydraulic root architecture models to macroscopic representations of root hydraulics (Vanderborght et al., 2021). The resulting model has basically three macroscopic parameters defined at the soil element scale, or at the plant scale, rather than for each segment of the root system architecture and are more appropriate for comparison with geophysical segments.
The root parameterization employed in hydrological models defines the "state model" which is commonly described by a state vector x SW(t) (e.g., soil water content) with a Prior PDF function p(x SW ). This prior PDF quantifies what we know about the state of the system given the prior distribution of the parameters and before considering the observed data.
Equation 1 describes the state vector at time t, F( . ) is the evolution operator, m is a model parameter vector with a Prior distribution p(m), and ω t is an uncertainty model error (Linde, 2014;Ghorbanidehno et al., 2020). For instance, soil moisture dynamics are described by the Richards' equation, which is implicitly contained into F( . ) and is highly non-linear. The model uncertainty ω t accounts for structural model errors (e.g., misrepresentation of physical processes), errors due to the discretization, parameter errors, initial condition errors, etc.
In general, the actual value of each error is unknown, thus it is often considered as a random variable described by a Gaussian distribution. One more pending issue in the application of RWU models is linked to the uncertainty and difficulty of measuring or predicting Root Length Density (RLD) with depth over time (Cai et al., 2018). Large uncertainties on the RLD definition are likely to lead to wrong estimates of model states. We discuss further how to reduce this uncertainty as a perspective in section Assimilation of Stem-Based Method Observation Data. What about the microscopic models? There isn't a lot mentioned about them.

Formulating a Soil/Root Electrical Conductivity Model and the Petrophysical Relationship
Electrical resistivity tomography produces a set of injected electrical currents and corresponding voltages. Poisson's equation, which expresses the relationship between the measured voltage and the electrical conductivity of the medium, given the injected current, is used to forward model the electrical potential data (Binley and Kemna, 2005). Electrical conductivity is a property of a medium that quantifies how easily electrical charges are conducted and is linked to the soil properties according to specific petrophysical relationships. The most widely used is Archie's law (Archie, 1942), which empirically defines the dependence of the bulk soil electrical conductivity on porosity, pore water conductivity and soil saturation). Nevertheless, petrophysical relationships are site specific and need to be calibrated before serving as translation between ER and SMC. This implies that petrophysical relationships, such as Archie's law, have some parameters to be calibrated. A review of the existing methods and limitations is presented by Cimpoiaşu et al. (2020). Michot et al. (2016) showed that the petrophysical relationships in the root zone are biased due to the influence of the root system on the measured ER. In a synthetic study, Rao et al. (2019) evaluated how geometrical anisotropy of root networks impact the effective conductivity of the rhizosphere (the zone near the root). As petrophysical model uncertainty can be influenced by the presence of a root system, a better knowledge of root electrical properties would help define the formulation of a mixed root soil ER in order to refine the petrophysical relationships (Ehosioke et al., 2020).
The impact of roots on the impedance properties of the rhizosphere has been observed and modeled, even for relatively small fractional volumes (Rao et al., 2019). This impact can be direct and indirect. The former derives from the direct rootcurrent iteration and depends on the large surface area and capacitance of roots, which, in turns, has been linked to root physiology and suberization (Weigand, 2017;Peruzzo et al., 2021;Tsukanov and Schwartz, 2021). The latter, indirect impact, is a combination of the root-induced changes in the physical, chemical, and biological properties of the rhizosphere, e.g., RWU, root exudation, solute uptake and rhizosphere microbiome changes (Landl et al., 2021).
The translation (or mapping) of soil and root electrical resistivity to SWC defines the so-called 'observation model'. The observation model (i.e., predicted ER observations) is commonly described by a state vector y ER(t) at assimilation time t with a prior PDF function p(y ER(t) ). H ER ( . ) denotes the observation operator describing the Poisson's equation, x ER(t) defines the ER values after petrophysical transformation of x SW(t) , n is a parameter vector, and v t is a stochastic model error.
The term v t incorporates both measurement errors and observation model uncertainties, by taking into account the following criteria: -Petrophysical model uncertainty should be considered in order to avoid the creation of unrepresentative hydrological models from geophysical inversions (Linde, 2014;Brunetti and Linde, 2018;Tso et al., 2019). The Archie's law parameters can be assumed as stochastic multivariate Gaussian distribution vectors (Botto et al., 2018) or as multivariate uniform distributions (Tso et al., 2020). -Regarding the ERT data, Rossi et al. (2015) stated that the likelihood functions can be obtained from the measurement error PDFs defined as a zero-mean Gaussian measurement with a given variance. -The uncertainty related to the mismatch between the spatial resolution of the geophysical parameters and the scale that characterizes the hydrogeological parameters should be considered.

DEGREE OF INTEGRATION BETWEEN DATA AND MODEL
The degree of integration of data with models may range from simple trial-and-error calibration (e.g., Binley et al., 2002) to fully coupled hydrogeophysical inversion and DA (e.g., Rajabi et al., 2018). Huisman et al. (2012) discussed whether the information content of the data is sufficient to obtain reliable estimates of model parameters. However, in general, the availability of spatially extensive (and time intensive) data greatly improves the model capability to identify the relevant governing hydraulic parameters, and/or unknown forcing conditions. There are several methods to compare quantitatively observed and simulated data such as moment analysis (e.g., Binley et al., 2002;Monego et al., 2010;Haaken et al., 2017) or the datadomain correlation approach of Johnson et al. (2009) that focuses on maximizing the correlation of hydrological and geophysical time series. Binley et al. (1996) introduced the concept of pixel breakthrough curves for analyzing time-lapse geophysical inversions of solute transport experiments. One may use such system responses to calibrate a hydrological model (e.g., Kemna et al., 2002;Crestani et al., 2015). For instance, Mary et al. (2020) used ERT and Mise-à-la-masse (MALM) data in order to constrain the initial values and boundary conditions of hydrological models. Cassiani et al. (2016) applied hydrological modeling for a quantitative interpretation of the ERT experiment on a plant-root system. The next degree of integration consists in incorporating additional knowledge during two sequential inversions, i.e., applied to the ERT measurements and the Richards' equation. This latter approach does not ensure an accurate quantitative description of the physical state, typically violating the mass balance  as a result of resolution limitation in geophysical data. Conversely, a coupled inversion links a hydrological model directly with a forward model of the geophysical data, and the mismatch between measured and modeled geophysical response is minimized (Wagner et al., 2019;Wagner and Uhlemann, 2021). In doing so, the soil hydraulic parameters used in the hydrological model can be optimized, while error propagation is avoided (Rubin and Hubbard, 2005;Vereecken et al., 2006;Binley et al., 2010Binley et al., , 2015Hinnell et al., 2010). Geolectrical inversion has known issues, typically due to the non-uniqueness of the solution in addition to incomplete or imperfect data due to practical limitations in collecting them. Advantages and drawbacks (structural errors in the hydrological conceptual model, dependency on known petrophysical relationships) of coupled inversions have been discussed by many authors (e.g., Hinnell et al., 2010;Camporese et al., 2015a;Slater and Binley, 2021;Yu et al., 2021). An emerging method for combining model predictions with observations is provided by the DA framework. The biggest advantage of DA methods over classical inversion is that they do not require inverted ER values, thus reducing the uncertainty in the estimation of the hydrological quantities, since no artifacts are inserted in the method by solving a classical inverse problem (see section Model Parametrization). Several studies show the promises of DA schemes for dynamic model parametrization in the context of soil-plant modeling (e.g., Han et al., 2012;Hu et al., 2017Hu et al., , 2019.

MODEL PARAMETRIZATION Prior Models and Global Sensitivity Analysis
As a choice for the prior information we generally want to define minimally subjective Prior PDFs (see sections Root water uptake models, Formulating a Soil/Root Electrical Conductivity Model and the Petrophysical Relationship). A Gaussian prior will enforce smoothness via the prior covariance. Hence, Gaussian priors are not suitable for characterizing discontinuous (non-smooth) fields and we would prefer edge-preserving priors (Arridge et al., 2019). Note that prior information is also useful to constrain the deterministic inversion by setting the lower and upper bounds for each parameter (e.g., Kuhl et al., 2021). Prior knowledge can be estimated from onsite calibration measurement or from literature values. The number of parameters governing a hydrological system can be high. A general rule of thumb is that parameters that are not being updated are assumed to be known, and literature values are used throughout the DA process. In order to decide which of the model and observation parameters need to be fixed, a global sensitivity analysis is often useful. To this end, Kuhl et al. (2021) quantified the sensitivity of each tested value by computing the root mean square error (RMSE) between the reference synthetic ER data and the perturbed ER data modeled. More advanced approaches can be used such as the Morris' method (1991) employed to determine the few most influential parameters among a large number of parameters .

Definition of the Objective Function
Both for deterministic coupled hydrogeophysical inversion and DA methods, an objective function (or cost function) is defined, which expresses the misfit between observations and model predictions. This function can take different forms when considering the hydrogeophysical inversion (Rubin and Hubbard, 2005;Vereecken et al., 2006;Hinnell et al., 2010;Mboh et al., 2012).
To express the objective function for the variational approach of DA, we follow the nomenclature used by Ning et al. (2014). In the general case, we aim at estimating a time trajectory of the state vector x = [x 0 , x 1 . . . , x T ] (with the time steps ranging from 1 to T and T being the assimilation window), by minimizing the cost function: The cost function is explicitly expressed as the sum of three components to draw the similarities with classical inversion of hydrogeophysical data. J b (x), J q (x) and J 0 (x) stand respectively for the a priori, the model and the data term. Here, x b is a prior estimate of the true initial state commonly referred to as the background state. The notation e 2 P −1 b represents the weighted quadratic norm, that is, e T P −1 e, where T is the transpose operator, and P represents the covariances of the error e (same apply for e 2 Q −1 t and e 2 R −1 t where Q and R represent the covariances of the error e for the model and the data weighted quadratic norm, respectively). In section Optimal Mass Frontiers in Water | www.frontiersin.org Transport Metric in Variational Data Assimilation we introduce a new objective function based on a different norm in order to improve the estimation of the misfit for RWU models.

Model Update
After the prediction step (Eqs. 1, 2), DA requires one additional step called update. This sequence is repeated recursively to estimate unknown states at multiple times of the system's evolution when observations are available. Note that this constitutes the main difference with single-step (also termed "batch") geostatistical inversion where all the time steps are assimilated at once. Among DA approaches, we distinguish between Bayesian filters and Bayesian smoothers. While Bayesian filters (e.g., Kalman Filter, KF) in its basic form only computes estimates of the current state of the system given the history of measurements, Bayesian smoothers (e.g., Kalman Smoother, KS) can be used to reconstruct states that happened before the current time. Both filters and smoothers suffer from the computational cost related to the propagation over time of the covariance matrices. To overcome this limitation, the covariance matrices can be reduced by representing them with a number of random realizations generated from an appropriate distribution (Ensemble based methods). The Kalman filter can be extended to the Ensemble Kalman filter (EnKF) which allows for nonlinear models but is limited by the Gaussian assumption. An alternative is the Particle Filter (PF) approach which extends the application of DA methods to non-Gaussian parameters. Many references report successful applications of both smoother and filter families to integrate geophysical observations (Pasetto et al., 2012;Manoli et al., 2015;Rossi et al., 2015) but also plant and atmospheric observations (Dong et al., 2016) into hydrological models.

TOWARD AN INTEGRATED RWU MODEL USING DA
DA algorithms rely on the state model, the observations and the observation models (petrophysical relationships) described already in sections Root Water Uptake Models and Formulating a Soil/Root Electrical Conductivity Model and the Petrophysical Relationship. For the specific case of improving the combination of geoelectrical data with hydrological root model we propose two pathways for improvement.

Optimal Mass Transport Metric in Variational Data Assimilation
In section Definition of the Objective Function, we have seen that in most cases the discrepancy between measured and predicted data is expressed as the Euclidean L 2 misfit function. Nevertheless, in reality, as the root system forms a network this norm is not optimal to capture its complex geometry. Recent studies show the potential of the optimal transport theory in defining a more suitable objective function based on the Wasserstein distance (Eq. 4, based on the formalism of Ning et al., 2014). Several studies report the advantages of using Wasserstein distance compared to Euclidean distance. For example, Ning et al. (2014) advocate the relevance of transportation metrics for quantifying non-random model error in variational DA for non-negative natural states and fluxes when Prior knowledge is not always reliably defined. Similarly, Métivier et al. (2016) use the Kantorovich-Rubinstein (KR) norm as an optimal transport distance to improve inversion of full waveform seismic tomography while Hernández and Liang (2018) propose a hybrid Bayesian and variational data assimilation method for highresolution hydrologic forecasting.

Assimilation of Stem-Based Method Observation Data
Most of the studies cited above consider empirical functions describing the root system length or density needed for the macroscopic RWU model. Stem based geoelectrical methods such as electrical capacitance or MALM methods aim at retrieving the RLD via the distribution of current source density (CSD) (Dalton, 1995;Peruzzo et al., 2020). In order to reduce the solution's non-uniqueness we can use: (i) a simplified observation model approach which consists in informing the root-zone sink term used in macroscopic RWU model with the known RLD inferred from inverted CSD data (see e.g., Mary et al., 2019), (ii) a more accurate model approach with a forward model predicting the CSD based on the solution of the Poisson equation. This DA approach would allow an assimilation of either ERT or both ERT and CSD datasets, taking into account the associated measurement uncertainties.

CONCLUSIONS
Time lapse ERT measurements are the most informative tools to derive soil water movement and variation of soil water content due to root water uptake. Thus, it appears natural to cast geoelectrical monitoring data within a dynamical parametrization framework. This is especially straightforward considering the non-linearity and complexity of the soil-plant system and all the uncertainties that need to be accounted for. Improving our understanding of soil-plant interactions requires the development of coupled numerical hydrogeophysical models and data assimilation methods. Indeed, the governing parameters of soil-plant systems are likely non-Gaussian variables and the forward and translation models are highly non-linear. The parametrization strategy must then be adapted from coupled optimizers, specific DA schemes (Particle Filters) and/or improved error estimations (Wasserstein). Additionally, the recent advances in characterizing roots using geoelectrical methods such as stem based ones can be seen as a leverage to introduce more extensive root and soil knowledge into the hydrological model. We advocated the use of macroscopic RWU models which have the advantage of being easily informed by recent development in root imaging (CSD). Considering that one of the biggest challenges is to evaluate uncertainties in the model parameters, geophysical data and forcing data, we can conclude that there is an obvious advantage of the Bayesian approach over deterministic methods.

AUTHOR CONTRIBUTIONS
BM: writing-original draft. All authors contributed equally to the manuscript writing-review and editing.

ACKNOWLEDGMENTS
We wish to acknowledge support from the project ECZ-Dry: New technologies to monitor the Earth Critical Zone in water-limited ecosystems funded by Italy-Israel Scientific and Technological Cooperation Programme (Scientific Track 2018). BM acknowledges the financial support from European Union's Horizon 2020 research and innovation programme under a Marie Sklodowska-Curie grant agreement (Grant no. 842922).