Time-Lapse ERT, moment analysis, and numerical modeling for estimating the hydraulic conductivity of unsaturated rock

Citation: De Carlo, L.; Farzamian, M.; Turturro, A.C.; Caputo, M.C. Time-Lapse ERT, Moment Analysis, and Numerical Modeling for Estimating the Hydraulic Conductivity of Unsaturated Rock. Water 2023,15, 332. https://doi.org/ 10.3390/w15020332 Academic Editor: M. Levent Kavvas Received: 29 November 2022 Revised: 29 December 2022 Accepted: 6 January 2023 Published: 12 January 2023 Copyright: © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). water Article Time-Lapse ERT, Moment Analysis, and Numerical Modeling for Estimating the Hydraulic Conductivity of Unsaturated Rock Lorenzo De Carlo 1,* , Mohammad Farzamian 2,3 , Antonietta Celeste Turturro 1 and Maria Clementina Caputo 1 1Water Research Institute, National Research Council of Italy, 70132 Bari, Italy 2Instituto Nacional de Investigação Agrária e Veterinária, 2780-157 Oeiras, Portugal 3Centre of Geographical Studies (CEG), IGOT, Universidade de Lisboa, 1600-276 Lisbon, Portugal *Correspondence: [email protected] Abstract: In recent years, geophysical techniques have been increasingly used to monitor flow and transport processes in the Earth critical zone (ECZ). Among these, electrical resistivity tomography (ERT) is a powerful tool used to predict hydrological parameters and state variables that influence the mentioned processes in the vadose zone because of the strong correlation between electrical and hydrological properties of the filtering medium. There have been many field tests considering geophysical prospecting in soils, where point scale hydrological sensors measurements are typically collected through sensors for geophysical data validation; on the contrary, when the unsaturated zone is made of hard rocks, the installation of such sensors is not a trivial issue owing to the extreme difficulties to guarantee contact between sensors and the surrounding medium. In this context, the geophysical data combined with appropriate numerical analysis techniques can effectively overcome the lack of information of the unsaturated subsurface, which is otherwise unpredictable with traditional methods. In the proposed case study, hydrogeophysical data were collected to provide a quantitative estimation of the hydraulic conductivity of sandstone through an integrated approach based on the moment analysis technique and numerical modeling. Keywords: unsaturated zone; hydrogeophysical monitoring; hydraulic conductivity; moment analysis 1. Introduction Hydraulic conductivity is the key parameter for the knowledge and understanding of the flow and transport processes that take place in the unsaturated zone. Many processes that involve groundwater management and protection, i.e., the evaluation of the aquifer recharge rate or the assessment of aquifer vulnerability from toxic substances, are strictly related to the hydraulic conductivity of the filtering medium. At the same time, the agricultural management practices depend on the accurate knowledge of the hydraulic properties of the soil. As known, the hydraulic conductivity of an unsaturated medium is a function of the water content or the water potential [ 1 ]. Few direct measurements of hydraulic conductivity function are carried out in a laboratory. The steady-state method was performed using a permeameter [ 2 ], while the hydraulic conductivity function was determined using the transient method [3]. In field conditions, natural processes such as ponded infiltration, copious rainfall, or irrigation lead to a near-saturated state owing to the air trapped in a significant fraction of the pore space, involving the field-saturated hydraulic conductivity, Kfs [4]. Infiltrometer tests were tested on both soils and rocks to provide a reliable estimation of such a parameter. Different types of infiltrometers were used on soil to measure K fs [ 5 – 8 ], while several criticisms increase the difficulties in performing such tests on rocks [ 9 – 11 ] because of the gap effect between the fragile probes and the rigid medium that could affect the uncertainty of the measurements. Water 2023,15, 332. https://doi.org/10.3390/w15020332 https://www.mdpi.com/journal/water Water 2023,15, 332 2 of 13 However, such experimental approaches do not take into account the wide spatial variability of the hydrological properties, assuming small-scale areas as representative of the investigated medium, as well as disturbances and errors that can affect the measurements. In some cases, unconventional large infiltrometer apparatus have been used to include as far as possible the heterogeneities of the subsurface [ 12 ] but, as a matter of fact, the problem of hydraulic characterization over large area still remains a great challenge for the scientific community. In the last decade, non-invasive geophysical techniques are increasingly used for hydrological issues because they have proven to be reliable, as well as relatively quick, inexpensive, and capable of investigating large areas with high data density sampling. Among geophysical techniques, electrical resistivity tomography (ERT) allows for imaging static electrical structures in the subsurface through resistance measurements collected on the ground surface or boreholes. A wide range of applications concern the use of ERT in the environmental [13], hydrological [14–16], and agricultural context [17,18]. In recent years, ERT surveys performed in time-lapse mode, i.e., repeating resistivity data collection along the same profiles or area at different temporal steps, have been revealed as a powerful tool for monitoring hydrological processes. The basic conceptual analogy between hydrological and electrical models led to the implementation of hydro-geophysical approaches, typically uncoupled and/or coupled, in order to constrain the hydrological model with geophysical data and to provide a quantitative estimation of the hydraulic conductivity and state variables [ 19 – 26 ]. The ERT-derived hydrological state variables are used to calibrate and constrain a hydrological model, thus identifying the corresponding governing parameters (i.e., hydraulic conductivity) with a wide range of numerical approaches, both stochastic and/or deterministic. Prediction of flow and solute transport was defined in terms of probabilistic properties of concentration [ 27 ]. In hydrology, several case studies considered the moment analysis approach to provide a measure of plume condition, as well as to evaluate forecasting scenarios in a long-term monitoring program [ 28 , 29 ]. In environmental features, geophysical datasets were increasingly used as input parameters for stochastic moment analysis [ 30 , 31 ], owing to the high capability to integrate hydrological data in areas not investigated with traditional measurements. In the proposed case study, a hydrogeophysical dataset, based on falling head infiltration measurements combined with a time-lapse ERT survey, was processed through a combined stochastic and deterministic approach in order to provide an initial accurate estimation of the saturated hydraulic conductivity, K s , of rock in unsaturated conditions. Quantitative information about the water mass and the movement of the center mass of the water plume were extracted from ERT-derived water content observations by performing the moment analysis technique. At the same time, water content distribution in the unsaturated zone during the infiltration test has been predicted through several simulation runs of the falling head test by constraining the water head and varying the K s value. The moment analysis was performed for each simulation and compared with the ERT-derived calculations, in order to provide an accurate estimation of K s . The goals of this paper are as follows: (1) to evaluate the ability of the geophysical tool to gain quantitative information of the dynamics occurring in the unsaturated rock at field scale, otherwise unpredictable with traditional observations; (2) to verify the reliability of stochastic techniques, such as moment analysis, in the monitoring of unsaturated flow; and (3) to reduce the uncertainty of such predictions by integrating geophysical datasets into stochastic and deterministic approaches. 2. Materials and Methods The approach used in the present paper combines a stochastic and deterministic processing of a hydrogeophysical dataset for estimating the K s parameter in the unsaturated zone. On one hand, the time-lapse ERT datasets collected during the infiltration test were used to image the resistivity variations of the subsurface that, in turn, were converted Water 2023,15, 332 3 of 13 into water content through Archie’s parameters calibrated in the laboratory. On the other hand, the unsaturated flow of the infiltration dynamics was simulated to predict the water content distribution in the rocky subsurface. The moment analysis calculated on both the ERT-derived water content dataset and the simulated one was used for connecting the two different approaches in order to provide an estimation of the K s . The conceptual scheme of the proposed approach is described in the flowchart shown in Figure 1. Water 2023, 15, x FOR PEER REVIEW 3 of 14 2. Materials and Methods The approach used in the present paper combines a stochastic and deterministic processing of a hydrogeophysical dataset for estimating the K s parameter in the unsaturated zone. On one hand, the time-lapse ERT datasets collected during the infiltration test were used to image the resistivity variations of the subsurface that, in turn, were converted into water content through Archie’s parameters calibrated in the laboratory. On the other hand, the unsaturated flow of the infiltration dynamics was simulated to predict the water content distribution in the rocky subsurface. The moment analysis calculated on both the ERT-derived water content dataset and the simulated one was used for connecting the two different approaches in order to provide an estimation of the K s . The conceptual scheme of the proposed approach is described in the flowchart shown in Figure 1. Figure 1. Workflow diagram describing the integrated hydrogeophysical approach used in the present case study for estimating saturated hydraulic conductivity, K s . 2.1. Properties of the Investigated Rock The rock considered in the case study is a Pleistocenic calcarenite, a sedimentary carbonatic porous sandstone of marine origin, widespread in the Mediterranean unsaturated zone. In terms of hydraulic properties, the rock is well characterized, mainly in laboratory activities, as described by an extensive literature. Caputo and Nimmo [32] provided accurate estimations of water retention and hydraulic conductivity functions of calcarenite samples collected in different rock quarries. Turturro et al. [33] tested the Darcy–Buckingham law in the centrifugal field on porous Figure 1. Workflow diagram describing the integrated hydrogeophysical approach used in the present case study for estimating saturated hydraulic conductivity, Ks. 2.1. Properties of the Investigated Rock The rock considered in the case study is a Pleistocenic calcarenite, a sedimentary carbonatic porous sandstone of marine origin, widespread in the Mediterranean unsaturated zone. In terms of hydraulic properties, the rock is well characterized, mainly in laboratory activities, as described by an extensive literature. Caputo and Nimmo [ 32 ] provided accurate estimations of water retention and hydraulic conductivity functions of calcarenite samples collected in different rock quarries. Turturro et al. [ 33 ] tested the Darcy–Buckingham law in the centrifugal field on porous rock media by demonstrating experimentally the direct proportionality between flux density and different centrifugal forces maintaining constant the water content. Turturro et al. [ 34 ] derived the whole retention curves of calcarenite rocks by a mercury intrusion porosimetry test for most of the mercury content range, filling the gaps exhibited by the retention curves obtained using the standard procedures. Caputo et al. [ 35 ] tested several laboratory methods for the prediction of the WRCs parameters critical for the numerical simulation of the unsaturated flow. Andriani et al. [ 36 ] conducted falling head infiltration tests at bench Water 2023,15, 332 4 of 13 scale on regularly shaped calcarenitic rocks with the aim to define the hydraulic behavior of calcarenites with different textures in unsaturated conditions. 2.2. The Hydrogeophysical Test The experimental test considered in the present paper concerns a 4D time-lapse ERT monitoring of infiltration measurements performed at falling head on the bottom of a calcarenite quarry. Details of the experimental set up of the hydrogeophysical test, as well as elaborations of the infiltration test and geophysical results, were described in detail in [ 21 ]. The test consisted of a double injection of about 7 L of water in two different temporal steps into an infiltrometer ring of 0.50 m in diameter sealed to the outcropped rock. During the test, the decreasing of water level was measured with a pressure transductor (PTX DRUCK LTD, Stoke Prior, Bromsgrove, UK) and the deepening of the wetting front was monitored by 4D time-lapse ERT surveys. Twenty-five partial hydraulic heads, defined as recharge periods in the code, measured at specific times points during the infiltrometer test, were used as constrained boundary conditions for the simulations. Table 1reports the main information of the infiltration test, including time of the start and stop of the two injections, hydraulic heads, and ERT observation. For simplicity, only ten hydraulic heads are mentioned, corresponding to the values measured at the start and stop of the injections, and those observed at six ERT observations. Table 1. Time points of the infiltrometer test and ERT data collection. Infiltration Test Time Point (hh:mm) Hydraulic Head (cm) ERT Observation 11:20 t0 Start first injection 11:43 3.1 12:03 2.0 t1 12:33 1.0 t4 Stop infiltration measurements 12:43 0.8 13:23 t8 Start second injection 14:01 3.3 14:16 2.5 t12 14:36 2.1 t13 14:59 1.6 t14 Stop infiltration measurements 15:09 1.3 15:59 0.2 t15 In this paper, the geophysical outputs are embedded into a stochastic and deterministic scheme, as part of the integrated data processing described in Figure 1. The experimental set up of the hydrogeophysical test is shown in Figure 2. Forty-eight stainless steel electrodes were placed around the ring with 0.40 m interelectrodes spacing to collect about 1200 dipole–dipole resistivity data points with Syscal Pro Switch 48 (Iris Instruments) resistivity-meter. The data collection was repeated 15 times over 4.5 h in order to detect changes in resistivity caused by water content variations. Each acquisition time took about 15 min. ERTLab commercial software (Multi-Phase Technologies, LLC, in collaboration with GeoStudi Astier s.r.l.) was used for the inversion process. According to [ 25 ], a convenient formulation of Archie’s law was used to convert the inverted resistivity into saturation degree. ρ(ti)=ρ(t0)Sw(t0) Sw(ti)n , (1) where ρ (t 0 ) and ρ (t 1 ) are the resistivities of the rock ( Ω· m) at times t 0 and t 1 , respectively, and S w (t 0 ) and S w (t 1 ) are the saturation degrees at times t 0 and t 1 , respectively. This approach allows the simplification of Archie’s equation, as the “n” saturation index is Water 2023,15, 332 5 of 13 the only unknown parameter to be calibrated. Archie’s calibration was performed in a laboratory on calcarenite core samples. In the laboratory, starting from the saturation condition, the resistivity–water content curve was recorded every five minutes during the drying process. Water 2023, 15, x FOR PEER REVIEW 5 of 14 process. According to [25], a convenient formulation of Archie’s law was used to convert the inverted resistivity into saturation degree. ρ(t)=ρ(t)󰇧S(t) S(t)󰇨, (1) where ρ(t0) and ρ(t1) are the resistivities of the rock (Ω∙m) at times t0 and t1, respectively, and Sw (t0) and Sw (t1) are the saturation degrees at times t0 and t1, respectively. This approach allows the simplification of Archie’s equation, as the “n” saturation index is the only unknown parameter to be calibrated. Archie’s calibration was performed in a laboratory on calcarenite core samples. In the laboratory, starting from the saturation condition, the resistivity–water content curve was recorded every five minutes during the drying process. Figure 2. Experimental set up of the infiltrometer test: (a) technical drawing of the experimental set up in a plan view; (b) picture of 3D electrodes’ configuration arrangement. 2.3. Forward Hydrological Modeling: Richards’ Equation VS2DRTI code [37,38] was used for simulating the falling head infiltration experiment. The code solves the two-dimensional approximation of Richards’ equation over time using the finite-difference method, in both cartesian and radial coordinates, thus providing a quasi-three-dimensional model of an isotropic system, assuming angular symmetry. The 3D Richards’ equation for variably saturated flow in a homogenous rigid porous medium [2] is expressed in the following form: 𝜕𝜃(ℎ) 𝜕(𝑡)=𝜕 𝜕𝑥𝐾(ℎ)𝜕ℎ 𝜕𝑥+𝜕 𝜕𝑦𝐾(ℎ)𝜕ℎ 𝜕𝑦+𝜕 𝜕𝑧𝐾(ℎ)𝜕ℎ 𝜕𝑧+1−𝑆, (2) where z is the elevation, t is time, θ is the volumetric water content, h is the water potential, θ(h) is the water retention function, K(h) is the hydraulic conductivity function, and S is a sink term representing root water uptake. According to the van Genuchten model [39], the soil hydraulic properties were described as follows: 𝜃=𝜃+𝜃−𝜃 󰇟1+|𝛼ℎ|󰇠, (3) where θ is the water content (m3∙m−3); θr is the residual water content (m3∙m−3); θs is the saturated water content (m3∙m−3); h is water potential (kPa); α is a scale parameter inversely proportional to mean pore diameter (cm−1); and n and m are the shape parameters of soil water characteristic, m = 1 − 1/n, 0 < m < 1. According to [35], van Genuchten’s parameters were set as reported in Table 2. Figure 2. Experimental set up of the infiltrometer test: ( a ) technical drawing of the experimental set up in a plan view; (b) picture of 3D electrodes’ configuration arrangement. 2.3. Forward Hydrological Modeling: Richards’ Equation VS2DRTI code [ 37 , 38 ] was used for simulating the falling head infiltration experiment. The code solves the two-dimensional approximation of Richards’ equation over time using the finite-difference method, in both cartesian and radial coordinates, thus providing a quasi-three-dimensional model of an isotropic system, assuming angular symmetry. The 3D Richards’ equation for variably saturated flow in a homogenous rigid porous medium [2] is expressed in the following form: ∂θ(h) ∂(t)=∂ ∂xK(h)∂h ∂x+∂ ∂yK(h)∂h ∂y+∂ ∂zK(h)∂h ∂z+1−S, (2) where zis the elevation, tis time, θ is the volumetric water content, his the water potential, θ (h) is the water retention function, K(h) is the hydraulic conductivity function, and Sis a sink term representing root water uptake. According to the van Genuchten model [ 39 ], the soil hydraulic properties were described as follows: θ=θr+θs−θr 1+|αh|nm, (3) where θ is the water content (m 3· m −3 ); θr is the residual water content (m 3· m −3 ); θs is the saturated water content (m 3· m −3 ); his water potential (kPa); α is a scale parameter inversely proportional to mean pore diameter (cm −1 ); and nand mare the shape parameters of soil water characteristic, m= 1 − 1/n, 0 < m< 1. According to [ 35 ], van Genuchten’s parameters were set as reported in Table 2. Other parameters are required to be set in the VS2DRTI code, such as the anisotropy ratio K hh /K zz (dimensionless), the specific storage S s (m −1 ), the effective porosity ϕ (m 3· m −3 ), and the initial moisture θcontent (m3·m−3). For the present study, K hh /K zz and S s were set equal to 1 and 1.6 × 10 −4 m −1 , respectively, although small variations do not cause significant changes in the model output. Moreover, effective porosity ϕ was set equal to 0.45 on the basis of previous tests performed on core samples and the initial content θ to 0.22 m 3· m −3 , as estimated from the ERT-derived value from Archie’s conversion. Water 2023,15, 332 6 of 13 Table 2. Parametrization of the model domain. Parameter Scenario A Scenario B Scenario C Scenario D Scenario E Saturated Khh (cm·min−1)0.1 0.25 0.35 0.75 1 Khh/Kzz 1 Specific storage, Ss(m−1) 1.6 ×10−4 Effective porosity, ϕ0.45 θ, initial moisture content (m3·m−3)0.22 θr, residual moisture content (m3·m−3)0.02547 α(cm−1)0.07721 n 1.7541 Therefore, once the van Genuchten’s parameters have been set in VS2DRTI for all simulation runs, the only unknown variable is the saturated hydraulic conductivity K s . Several simulations scenarios were run with K s ranging from 0.1 cm · min −1 to 1 cm · min −1 (Table 2), according to the expected values reported in the literature. A 2D flow domain of 7 × 1 m (length × depth) was built with a constant grid spacing in the vertical direction, equal to 0.02 m, and variable grid spacing in the horizontal direction, ranging from 0.062 m to 0.143 m below the infiltrometer ring. Thus, the domain was discretized into 2912 nodal points. A no-flow boundary condition (Dirichlet condition) was imposed along the lateral faces of the domain and top surface outside the infiltrometer ring, while a free drainage boundary condition was used at the bottom. Simulation output times correspond to the ERT observations, i.e., t 1 = 20 min, t 4 = 50 min, t 8 = 100 min, t 12 = 153 min, t 13 = 173 min, t 14 = 196 min, and t 15 = 256 min after the start of the first injection, as reported in Table 1. 2.4. Moment Analysis The basics of the moment analysis is the following equation: Mijk(t) = yΓ ∆θ(x,y,z,t)xiyjzkdxdydz (4) The zeroth and first spatial moments correspond to i+j+k= 0 and 1, respectively. ∆θ is the water content changes based on the resistivity changes estimation inferred from the time-lapse ERT model, after removing the background water content. Γ is the volume of interest. Assuming angular symmetry of the monitored infiltration dynamics, a 2D modeling was analyzed in this case study to simplify the calculation of the moment analysis. The 2D section considered for the moments calculation crosses the ring through its center and, for a precise comparison of the different datasets, it corresponds to the same cross section where the ERT-derived water content was estimated. Therefore, the zeroth moment, M 00 , is the changes in water mass within the domain respect to the background (Equation (5)) and represents the water storage along the reference section, expressed in m3m−3. M00(t) = x∆θ(x,z)dxdz. (5) The first moment, M 01 normalized by the mass M 00 , defines the vertical center of mass of the plume at a given time, z, expressed by Equation (6). z=M01 M00 (6) Water 2023,15, 332 7 of 13 3. Results 3.1. ERT-Derived Water Content Outputs Figure 3shows the ERT-derived water content distribution during the infiltration test inferred at several time points. Water 2023, 15, x FOR PEER REVIEW 7 of 14 𝑀(𝑡)=𝛥𝜃(𝑥,𝑧)𝑑𝑥𝑑𝑧. (5) The first moment, M01 normalized by the mass M00, defines the vertical center of mass of the plume at a given time, z, expressed by Equation (6). 𝑧=𝑀 𝑀 (6) 3. Results 3.1. ERT-Derived Water Content Outputs Figure 3 shows the ERT-derived water content distribution during the infiltration test inferred at several time points. The background water content (Figure 3a) highlights a narrow water content range in the subsurface, from 0.20 to 0.22 m3∙m−3 estimated from Archie’s conversion, by denoting an almost homogeneous initial condition of the subsurface, as expected. When the first injection starts, water infiltrates below the ring, deepening over time until the end of the first injection (Figure 3b–d). The water content observed in the background conditions can be attributed to copious precipitation some days before the test, leading to high values of water content in the upper portion of the subsurface soon after the starting of the first injection. As shown in Figure 3e–h, the second injection increases the subsurface water content in the deeper layers of the subsurface (green fill turns yellow and then red). Figure 3. ERT-derived water content at different time points: (a) before the starting of the first injection; (b) 20 min; (c) 50 min; (d) 100 min after the first injection; (e) 25 min; (f) 45 min; (g) 68 min; (h) 128 min after the second injection. 3.2. Moment Analysis Derived from the ERT Dataset The zeroth moment, calculated from the ERT-derived water content outputs inferred along the 2D reference cross section, is shown in Figure 4. Figure 3. ERT-derived water content at different time points: ( a ) before the starting of the first injection; ( b ) 20 min; ( c ) 50 min; ( d ) 100 min after the first injection; ( e ) 25 min; ( f ) 45 min; ( g ) 68 min; (h) 128 min after the second injection. The background water content (Figure 3a) highlights a narrow water content range in the subsurface, from 0.20 to 0.22 m 3· m −3 estimated from Archie’s conversion, by denoting an almost homogeneous initial condition of the subsurface, as expected. When the first injection starts, water infiltrates below the ring, deepening over time until the end of the first injection (Figure 3b–d). The water content observed in the background conditions can be attributed to copious precipitation some days before the test, leading to high values of water content in the upper portion of the subsurface soon after the starting of the first injection. As shown in Figure 3e–h, the second injection increases the subsurface water content in the deeper layers of the subsurface (green fill turns yellow and then red). 3.2. Moment Analysis Derived from the ERT Dataset The zeroth moment, calculated from the ERT-derived water content outputs inferred along the 2D reference cross section, is shown in Figure 4. Soon after the first injection, the increase in water content below the infiltrometric ring, which leads to an increase in water storage, can be observed. As expected, the water storage gradually increases over time owing to the double injection. No significant difference in the slope of the curve can be observed after the second injection, probably because of the small change in water content below the infiltrometer ring and the short temporal window between the two injections. Water 2023,15, 332 8 of 13 Water 2023, 15, x FOR PEER REVIEW 8 of 14 Soon after the first injection, the increase in water content below the infiltrometric ring, which leads to an increase in water storage, can be observed. As expected, the water storage gradually increases over time owing to the double injection. No significant difference in the slope of the curve can be observed after the second injection, probably because of the small change in water content below the infiltrometer ring and the short temporal window between the two injections. Figure 4. Water storage estimated through the zeroth moment M 00 . Figure 5 shows the prediction of the movement of the mass center calculated using the water content values inferred from the time-lapse ERT outputs. No information can be observed within the first 0.20 m, which is about the spatial resolution of the ERT images. After the first injection, the mass center of the water bulb moves downward, reaching a maximum depth of 0.44 m from ground surface at time point t 8 . After the second injection, the added infiltrated mass water causes a rise in the mass center, as clearly observed in the shape of the curve, reaching a value of 0.35 m at time t 15 . Figure 5. Depth of the mass center reached after two injections. Figure 4. Water storage estimated through the zeroth moment M00. Figure 5shows the prediction of the movement of the mass center calculated using the water content values inferred from the time-lapse ERT outputs. No information can be observed within the first 0.20 m, which is about the spatial resolution of the ERT images. After the first injection, the mass center of the water bulb moves downward, reaching a maximum depth of 0.44 m from ground surface at time point t 8 . After the second injection, the added infiltrated mass water causes a rise in the mass center, as clearly observed in the shape of the curve, reaching a value of 0.35 m at time t15. Water 2023, 15, x FOR PEER REVIEW 8 of 14 Soon after the first injection, the increase in water content below the infiltrometric ring, which leads to an increase in water storage, can be observed. As expected, the water storage gradually increases over time owing to the double injection. No significant difference in the slope of the curve can be observed after the second injection, probably because of the small change in water content below the infiltrometer ring and the short temporal window between the two injections. Figure 4. Water storage estimated through the zeroth moment M 00 . Figure 5 shows the prediction of the movement of the mass center calculated using the water content values inferred from the time-lapse ERT outputs. No information can be observed within the first 0.20 m, which is about the spatial resolution of the ERT images. After the first injection, the mass center of the water bulb moves downward, reaching a maximum depth of 0.44 m from ground surface at time point t 8 . After the second injection, the added infiltrated mass water causes a rise in the mass center, as clearly observed in the shape of the curve, reaching a value of 0.35 m at time t 15 . Figure 5. Depth of the mass center reached after two injections. Figure 5. Depth of the mass center reached after two injections. 3.3. Numerical Simulations Figure 6shows the comparison between the ERT-derived water content (Figure 6a) and three different scenarios selected among five simulation runs. Particularly, scenario A underestimates the K s (Figure 6b), scenario E overestimates the K s (Figure 6d), and scenario C approximates the true distribution of soil moisture (Figure 6c). In order to provide a quantitative estimation of K s , the moment analysis for each simulation scenario was calculated for all five scenarios. Water 2023,15, 332 9 of 13 The comparison between the depth of the mass center estimated from the geophysical and numerical approaches shows that the ERT-derived K s value fits well the simulated one in the range 0.25 < Ks< 0.35 cm·min−1(Figure 7). Water 2023, 15, x FOR PEER REVIEW 9 of 14 3.3. Numerical Simulations Figure 6 shows the comparison between the ERT-derived water content (Figure 6a) and three different scenarios selected among five simulation runs. Particularly, scenario A underestimates the Ks (Figure 6b), scenario E overestimates the Ks (Figure 6d), and scenario C approximates the true distribution of soil moisture (Figure 6c). In order to provide a quantitative estimation of Ks, the moment analysis for each simulation scenario was calculated for all five scenarios. Figure 6. Comparison of the ERT-derived water content (a) and the simulated one for three different scenarios: (b) scenario A; (c) scenario C; and (d) scenario E. The comparison between the depth of the mass center estimated from the geophysical and numerical approaches shows that the ERT-derived Ks value fits well the simulated one in the range 0.25 < Ks < 0.35 cm∙min−1 (Figure 7). Figure 6. Comparison of the ERT-derived water content ( a ) and the simulated one for three different scenarios: (b) scenario A; (c) scenario C; and (d) scenario E. Water 2023, 15, x FOR PEER REVIEW 10 of 14 Figure 7. Comparison between the depth of the mass center estimated from ERT and simulations. 4. Discussion and Conclusions This paper presented an integrated approach based on hydrogeophysical data embedded into a stochastic and deterministic model scheme in order to provide quantitative hydraulic information on an unsaturated porous rocky subsurface. The understanding of the flow and transport processes in the unsaturated zone is an exciting scientific challenge and the prediction of hydraulic properties from geophysical data is an open question still to be fully solved. The added value pointed out in this manuscript is based on the capability of the integrated approach to infer unsaturated hydrological properties in a complex domain, when data are lacking, or where data collection is extremely difficult or affected by disturbances and errors. In this paper, the achievement of the objectives is obtained through a multi-step elaboration process. The first step is based on the geophysical data processing, which revealed the great potential of ERT in the characterization and monitoring of the unsaturated flow processes in the rocky subsurface, as shown in Figure 3. The ERT-derived water content cross sections showed a fast infiltration soon after the first injection (Figure 3b–d), followed by a slow increase in water content after the second injection (Figure 3e–h), owing to the high saturation degree of the rocky subsurface. Geophysical measurements alone can lead to misinterpretation, artifacts, or multiple models, being based on the intrinsic non-uniqueness of the geophysical solution. Therefore, in order to produce a robust hydrogeophysical model, in the second step, the geophysical outcomes were treated with a stochastic approach, based on the moment analysis technique, with the aim to track the movement of the injected water mass. In particular, the depth of the mass center was monitored during the infiltration dynamics (Figure 5). The geophysical tool is capable of detecting variations in the depth of the mass center: soon after the first injection, it increases up to about 0.40 m below ground surface (bgs), then it moves slightly upward after the second injection, reaching a depth of 0.35 m bgs at the end of the test. Changing of the mass center depth over time pointed out a detailed response of the unsaturated domain, which cannot be obtained using other traditional methods. This is a 0 10 20 30 40 50 60 0 50 100 150 200 250 300 350 400 Depth from ground surface (cm) Time since start of injection (min) Ks=0.35 cm/min ERT derived Ks=0.1 cm/min Ks=0.75 cm/min Ks=0.25 cm/min Ks=1 cm/min Figure 7. Comparison between the depth of the mass center estimated from ERT and simulations.