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.

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 root-current 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_tincorporates both measurement errors and observation model uncertainties, by taking into account the following criteria:

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 (Hu et al., 2017).

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₀, x₁…, 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⁰(x) stand respectively for the a priori, the model and the data term. Here, x_bis a prior estimate of the true initial state commonly referred to as the background state. The notation ||e||Pb-12 represents the weighted quadratic norm, that is, e^TP⁻¹e, where ^T is the transpose operator, and P represents the covariances of the error e (same apply for ||e||Qt-12and ||e||Rt-12 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 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.

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 non-linear 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.

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² 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).

Equation 4 describes the minimization of the time trajectory of the state vector x = [x₀, x₁…, x_T]. Compared to Eq. 3, only the term corresponding to the measurement error is still quantified by a quadratic norm, under the assumption that observational errors (||e||Rt-12) are well-represented by additive (Gaussian) random noise. The other terms are based on the Wasserstein metric W.

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 high-resolution hydrologic forecasting.

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.