Inverse estimation of hydraulic parameters of soils with rock fragments

.


Introduction
Soils with rock fragments (RF > 2 mm) impose additional challenges in assessing their hydraulic properties and hydrological processes compared with soils without RF, mainly because of their high heterogeneity and the impossibility of collecting undisturbed soil samples. While worldwide research on several characteristics of these soils has increased in the last decades, a consensual understanding of the effects of RF on hydrological processes still demands further investigation (Zhang et al., 2016). One of the biggest challenges addressed in the last years is determining soil hydraulic parameters needed to model water flow and contaminant transport through the vadose zone of these soils. Direct laboratory measurements on undisturbed soil samples are complex and impractical depending on the amount of RF, which requires the use of altered porous media created by mixing RF with finer soil particles (Ma and Shao, 2008;Beibei et al., 2009;Ma et al., 2010;Parajuli et al., 2017;Naseri et al., 2019). This strategy allows mimicking the effects of RF on hydraulic functions, but replacing the natural porous system with an artificial one introduces bias in the soil hydraulic functions. Alternatively, an indirect determination of soil hydraulic parameters through inverse modeling has been increasingly applied over the last decades (e. g., Š imůnek et al., 1998;Lazarovitch et al., 2007;Jacques et al., 2012;Rashid et al., 2015;Pinheiro et al., 2019;Brunetti et al., 2019) and is considered a feasible technique to estimate hydraulic parameters from data collected during non-destructive in-situ experiments in soils with RF. However, challenges to field applications of inverse modeling to experiments in complex porous media like those of soils with RF (Weiler and Flühler, 2004) still remain a major research gap.
For any inverse modeling problem, optimized parameters depend on the suitability and quality of (i) the experimental design, which comprises the boundary conditions, locations and time resolution of measurement sensors, and the degree of accuracy of the experimental data; (ii) the suitability of the transient variably-saturated water flow model and hydraulic functions; and (iii) the robustness of the optimization algorithm as determined by its convergence to a global minimum (Hopmans et al., 2002). To prevent issues with numerical convergence, the number of parameters being simultaneously fitted can be reduced (Ritter et al., 2003;Abbasi et al., 2003). A suitable model option is the Hydrus-1D software (Šimůnek et al., 2013), which constitutes a robust, transient, variably-saturated flow model often applied on inverse modeling problems (e.g., Lazarovitch et al., 2007;Jacques et al., 2012;Rashid et al., 2015;Pinheiro et al., 2019;Brunetti et al., 2019). Thus, the main challenge, especially for soils with RF, is defining suitable soil hydraulic functions and performing well-designed experiments.
A suitable method for determining parameters of hydraulic functions for layered soils is the field (in-situ) drainage experiment, known as the internal drainage or instantaneous profile method (Hillel et al., 1972;Dikinya, 2005;Mavimbela and Van Rensburg, 2013). The method consists of quasi-saturating a bare soil in which sensors (probes to measure water contents and pressure heads) have been installed at selected depths. The soil surface is then covered with a plastic sheet to avoid any evaporation, and the soil is allowed to drain by gravity while the transient drainage process is monitored by simultaneous measurements of soil water contents and pressure heads. Efforts to apply inverse modeling in the field drainage experiment have focused on several issues as, for example, the definitions of the hydraulic head gradient at the bottom boundary (Romano 1993), the likelihood of non-uniqueness and instability of the inverse solution due to an increase in the number of optimized parameters for layered soils (e.g., Zijlstra and Dane, 1996;Ritter et al., 2003), and the combination of water content and pressure head measurements (Dikinya, 2005). All these issues demand further investigations, but evaluating the most suitable hydraulic functions for soil with RF is an important and still unresolved investigation gap.
The choice of the hydraulic function shapes could be guided by prior knowledge about the soil pore size distribution. However, almost no knowledge of pore size distribution is available for in-situ soils with RF. Multimodal (soil hydraulic functions, having more than one maximum in the soil water capacity curve and usually corresponding to a pore system composed of more than one kind of pores, e.g., structural and intra-aggregate pores) rather than unimodal soil hydraulic functions, with only one maximum in the soil water capacity curve, may be more suitable for soils with RF due to the presence of large pores and cracks around RF and small textural pores in the aggregated fine material between them (Durner, 1994;Ma and Shao, 2008;Ma et al., 2010;Zhang et al., 2016). Yet, multimodal functions increase the number of parameters to be optimized in the inverse problem, not only increasing the likelihood of non-uniqueness of the inverse solution, but also requiring measured water contents and pressure heads at a finer temporal resolution so that the flow model can capture the effects of different pore domains on hydraulic functions.
However, water contents and pressure heads in field drainage experiments without evapotranspiration vary only over the narrow wet range of soil hydraulic functions. The soil hydraulic parameters optimized considering only such data may thus not represent complete pore domains and adequately define the behavior of soil hydraulic functions over the dry range beyond the drainage phase. This limitation can be avoided and/or minimized by defining the inverse solution using additional data involving water retention in the dry range (Hopmans et al., 2018;Pinheiro et al., 2019). Such data can be measured in the laboratory, for example, using a dew point potentiometer, an instrument that is also suitable for soils with RF (Parajuli et al., 2017;Gubiani et al., 2021b). Thus, the shape of soil hydraulic functions can be better determined by combining both data sets (in the wet and dry range) in the optimization algorithm.
Avoiding ill-posed problems due to an increase in the number of soil hydraulic parameters may represent the most significant challenge when combining multimodal soil hydraulic functions with layered soil. For example, combining the van Genuchten bimodal water retention curve (Durner, 1994) with a soil profile with three soil horizons results in a vector of 21 parameters (24 when including Ks) to be optimized. A solution to this inverse problem with so many soil hydraulic parameters may become ill-posed (e.g., Ritter et al., 2003;Abbasi et al., 2003), mainly when local gradient-based optimization algorithms, such as the Levenberg-Marquardt scheme, are used (Šimůnek et al., 2012). To address this problem, some parameters may be independently estimated and their values fixed in the inverse solution, reducing the dimensionality of the inverse problem (Abbasi et al., 2003). Another strategy is to perform sequential optimization (Ritter et al., 2003;Abbasi et al., 2003), which repeatedly runs the optimization algorithm for selected parameters at each step.
Therefore, the main objective of this study is to determine the hydraulic parameters of soils with RF without disturbing the soil's structure while addressing the main investigation gaps discussed above. The additional objectives of this study are to determine suitable hydraulic functions for soils with RF and find a sequential optimization strategy to estimate their parameters using inverse modeling and data from field drainage experiments and dewpoint potentiometer measurements.

Instrumentation
Capacitance sensors with a factory calibration (ECH20 10HS Large Volume Soil Moisture Sensors by Meter, Inc.) to measure water contents were horizontally installed in the soil profiles at depths of 5, 20, 40, and 65 cm at Location 1, depths of 5, 20, 35, and 55 cm at Location 2, depths 5, 17, and 40 cm at Location 3, and 5, 20, 35, and 55 cm at Location 4. Additionally, soil water potential sensors with a factory calibration (Teros 21 by Meter Inc.) were installed at Locations 2 and 4 (at depths of 5 and 20 cm) and Location 3 (5 cm). Note that the sensors were not installed at the center of each layer (see Section 2.1), and some layers had more than one sensor. This is because the sensors' depths were defined according to other studies' needs on water balance. A detailed description of the use of water content and pressure head measurements is provided in Section 2.4.1. After the sensor installation on November 7, 2020, soil trenches were backfilled with the original soil materials of each layer, followed by compaction to minimize preferential flow (Fig. 2).

Drainage experiments
Four drainage experiments were performed at each location on November 30, 2020, December 3, 2020, April 25, 2021, and May 22, 2021. The area was planted with soybean until March 15 (during the first two experiments) and with oat after April 10 (during the last two experiments). During all experiments, plants were at their initial growth stage, and they did not interfere in the experimental determinations.
A circular area of approximately 2 m 2 (sensors at its center) was delimited with a 20-cm high plastic sheet, 10 cm above the soil surface and 10 cm buried. Water was ponded on the surface of the boarded soil, and an increase in the soil water content was monitored with 10 HS sensors every 2 min. Water supply was ceased after five similar readings of the highest water content were observed at all depths. When all ponded water had infiltrated, the soil surface was covered with a plastic sheet to prevent evaporation. Water contents and pressure heads were then measured over time while redistribution and drainage took place (Fig. 2).
The experiments were performed shortly (two or three days) after significant rainfall events when the initial water content was close to field capacity. This guaranteed that the wetted soil volume was surrounded by wet soil, avoiding considerable lateral flow during the drainage experiments. No rainfall occurred shortly after the onset of the experiments, allowing enough time for drainage rates to decrease to negligible values. The duration of observations was 50 h for the experiments on November 30, 2020 (a rainfall occurred after 50 h), and 180 h for the experiments on December 3, 2020, April 25, 2021, and May 22, 2021.

Definition of hydraulic functions and parameter optimization by inverse modeling
Parameter optimization was performed using the inverse routine of the Hydrus-1D model (Šimůnek et al., 2013). The one-dimensional Richards equation solved by Hydrus-1D for vertical drainage flow is given as follows: (1) where C is the water capacity (dθ/dh, cm − 1 ), θ is the volumetric water content (cm 3 cm − 3 ), h is the soil pressure head (cm), t is time (min), z is the vertical spatial coordinate (cm), and K is the hydraulic conductivity (cm min − 1 ). From several hydraulic functions available in Hydrus-1D to describe the θ-h-K relations in equation (1) we selected one of the most common, a unimodal water retention function (van Genuchten, 1980), and its parameters were optimized considering the measured and estimated θ and h in the wet range of the drainage experiment. Subsequently, we used an expansion of this unimodal function to a bimodal one (Durner, 1994) to reduce inconsistencies that would prevent the extrapolation of the unimodal function to drier domain. These strategies are described in detail in the next sections.

Parameter optimization of the unimodal water retention curve
The van Genuchten-Mualem functions (van Genuchten, 1980) were used to describe the θ-h-K relation in equation (1).
where θ s1 and θ r1 are the saturated and residual water contents (cm 3 cm − 3 ), α 1 (cm − 1 ) and n 1 (-) are fitting shape parameters, K and K s1 (cm min − 1 ) are the hydraulic conductivities for unsaturated and saturated conditions, S e1 is the effective saturation (-), and l (-) is the poreconnectivity parameter. Subscript "1" distinguishes the parameters optimized for a unimodal curve from those optimized for a bimodal curve. Hysteresis was disregarded in all simulations.
During the inverse modeling, the parameters θ s1 , θ r1 , α 1 , n 1 , and K s1 were sequentially optimized (see below), while the pore-connectivity parameter l was kept at 0.5 according to Mualem (1976). All soil profiles ( Fig. 1) were represented in Hydrus-1D by a discretized scheme with 101 equally spaced (1 cm) nodes. Each soil profile consisted of three materials to represent the sequence of soil layers (see Section 2.1). The upper boundary condition was set as a zero flux (a covered soil profile), and the lower boundary was set as free drainage, in agreement with field conditions where no water table or impeding soil layers were present within several meters of depth. A pressure head of − 10 cm was assumed as the initial condition along the entire domain. Although the Teros 21 sensor cannot accurately measure pressure heads higher than − 100 cm, a value close to saturation was considered a reasonable approximation of the field-saturated condition, close to the air-entry value for coarse-textured soils (Rawls et al., 1982). The suggested default values by the Hydrus-1D model were used for the water content error tolerance (0.001 cm 3 cm − 3 ), pressure head error tolerance (-1 cm), and the maximum number of iterations (1 0 0) for the numerical solution of the Richards equation. Hydrus-1D automatically selects an optimal time step, allowed to vary between the minimum and maximum limits set as 0.0144 min and 7200 min, respectively.
Parameters θ s1 , θ r1 , α 1 , n 1 , and K s1 were optimized by minimizing the objective function (Šimůnek et al., 1998): where q j *(x, t i ) represents the specific measurements of θ and h (measurement sets j) at time t i at location x, q j (x, t i , b) represents the corresponding model predictions of θ and h for the vector of optimized parameters b (θ s1 , θ r1 , α 1 , n 1 , and K s1 ), v j and w i,j are weights associated with a particular measurement set or point, respectively, m q is the number of different sets of measurements (2 in this case, for water contents and pressure heads), and n qj is the number of measurements within a particular measurement set. The weighting coefficients w i,j were all fixed at 1, giving the same weights to all data points in a particular measurement set. Furthermore, the weighting coefficients v j were fixed at (Šimůnek and Hopmans, 2002a,b): which implies that the objective function Φ (equation (4)) is the average weighted squared deviation between the measured and fitted values normalized by the measurement variances (σ 2 ). The Hydrus-1D model additionally reports various other statistics, such as Mean Weighted Error (ME), Mean Weighted Absolute Error (MAE), Root Mean Square Weighted Error (RMSWE), Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC). The Root Mean Square Weighted Error (RMSWE), Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC) are defined in Hydrus-1D as follows.
where the sum of squared weighted deviations Φ is calculated using (4), n is the number of data points, and m is the number of optimized parameters. The Akaike and Bayesian information criteria (Akaike, 1974(Akaike, , 1977 penalize for adding fitting parameters, i.e., the model with the smallest number of parameters (which minimizes AIC and BIC) is a preferred fitting model (Šimůnek and Hopmans, 2002b). Data for q j *(x,t i ) consisted of the time series of θ measured during the drainage experiments with the 10 HS sensors and h measured with the Teros 21 sensor at the end of the drainage experiments (locations of measurements are described in Section 2.2). Only the final measured h was used because Teros 21 readings are insensitive at wetter conditions. At the locations where no measurements of the pressure heads were performed, estimates of h were used. At Location 1, the pressure heads h at depths 5 and 20 cm were considered equal to the mean h value at corresponding depths at Locations 2 and 4. The pressure heads h at larger depths without measurements were considered the same or about − 3 to − 10 cm less than h at the overlying position. This strategy facilitated the setting of initial values for hydraulic parameters before running the automatic Hydrus-1D "Inverse solution" routine (see step 2 in the next section). Observations points were defined in the discretized profile scheme of Hydrus-1D to obtain estimated q j (x, t i , b) values corresponding to measured q j *(x, t i ) values.

Sequential optimization
A sequential calibration was implemented to reduce nonconvergence of the optimization algorithm and the likelihood of nonuniqueness when a large number of parameters had to be optimized (e.g., Hopmans et al., 2018;Ritter et al., 2003;Abbasi et al., 2003). The entire vector of 15 parameters (θ s1 , θ r1 , α 1 , n 1 , and Ks 1 for each of the three soil layers) was reduced to a few selected parameters to be calibrated in different steps as described below.
1. The initial estimate of θ s1 , θ r1 , α 1 , n 1 , and K s1 for each material (soil layer) was taken from the Hydrus-1D catalog for a loam soil (θ r1 = 0.078, θ s = 0.4, α 1 = 0.036 cm − 1 , n 1 = 1.56, K s1 = 1.7⋅10 -2 cm min − 1 ). 2. Before running the Hydrus-1D "Inverse solution" routine, most parameters were several times updated in Hydrus-1D by a trial-anderror procedure. Parameters were approximately refined by considering measured data of θ and visual inspection of the similarity between the measured and simulated time series of θ in the Hydrus-1D graphical interface. A good improvement in the fit between measured and simulated data was obtained by trial-and-error when considering θ s1 and θ r1 close to the initial and final values, respectively, of measured θ in the drainage experiments. Subsequently, α 1 and n 1 were also modified to improve the similarity between simulated and measured water content values. All these improvements were best achieved with K s1 fixed at 0.2 cm min − 1 . 3. Parameter optimization with the Hydrus-1D inverse solution routine was accomplished in sequential steps to optimize only a few optimized parameters in each step, thus minimizing ill-posed problems and avoiding non-uniqueness and unidentifiability. First, only α 1 and K s1 were optimized, and they were forced to have the same values for all soil layers. This was done by manually setting Kod = 2 (a code indicating that the same parameter value should be used for all soil layers) in the input file Fit.in. Then, the initial estimates of K s1 and α 1 were updated with their respective optimized values. Next, only θ s1 and θ r1 for all layers were optimized (with Kod = 1, allowing different parameter values in different layers). After that, the initial estimates of θ s1 and θ r1 for all layers were updated with optimized values. In the next step, α 1 , n 1 , and K s1 for all layers (with Kod = 1) were optimized and updated. Finally, a last minor refinement was done for θ s1 and θ r1 , optimizing them for all layers. Thus, the number of optimized parameters increased from two to six and then to nine, respectively, for the first, second, and third optimization steps. The number of optimized parameters in the third step (nine) was higher than usually recommended to avoid ill-posed problems (Ritter et al., 2003;Abbasi et al., 2003). However, it is expected that the initial trial-and-error steps followed by the first and second steps of inverse modeling located the search algorithm in the parameter space quite close to the global minimum of the objective function, making the third step of the inverse solution less prone to ill-posedness, even with nine parameters being optimized simultaneously. Restrictions were imposed on n 1 (1.1 ≤ n 1 ≤ 3.5) and θ s1 and θ r1 in cases of nonconvergence.

Expanding the unimodal to bimodal water retention curve
Although the simulated θ and h matched their respective measurements quite well (RMSWE < 0.05 in most cases), the optimization process made the unimodal water retention curve overestimate the soil residual water content. Optimized values of θ r1 ranged between 0.26 and 0.44 cm 3 cm − 3 . These values are inconsistent for the coarse-textured soils under investigation (Fig. 1). Overestimation in optimized values of θ r1 could be a consequence of using only the θ and h values measured in the wet range in the objective function.
We attempted to eliminate this problem by also considering the dry range θ and h data measured in the laboratory (h < − 5000 cm) with a dew point device (Gubiani et al., 2021b) in the objective function. This made optimized values of θ r1 decrease, but the agreement between simulated and measured θ and h for the drainage experiments became significantly worse. Therefore, keeping the unimodal θ-h-k relations with the parameters optimized using θ and h measured in the wet range would provide accurate drainage predictions, but its use to simulate fluxes in a drier soil (e.g., root water uptake, infiltration) would most probably lead to inaccurate results.
We, therefore, concluded that the Hydrus-1D algorithms for water flow and parameter optimization yielded unimodal water retention and hydraulic conductivity functions, which were unable to describe soil hydraulic properties for the entire saturation range. A likely reason is that the porous system of the soils evaluated in this study is more complex than the unimodal model would allow, which is possibly a consequence of the presence of coarse material (Fig. 1). Thus, the unimodal θ-h-k relations were expanded to bimodal functions (Durner 1994): where θ s1 , α 1 , n 1 , l, and S e1 are parameters described with equations (2) and (3), w 1 and w 2 (dimensionless) are weighing factors for each sub curve (conditioned to w 1 + w 2 = 1), and θ r2 , α 2 , n 2 , S e2 , and K s2 have the same meaning as their corresponding parameters in equation (2) and (3), but now for smaller textural pores. Optimizing all parameters of equations (9) and (10) with inverse modeling would likely be unfeasible and highly dependent on their initial guesses. We assumed that the best initial guess could be provided by keeping θ s1 , α 1 , n 1 , and m 1 at their previously optimized values (Section 2.4.2) and estimating w 1 , w 2 , θ r2 , α 2 , and n 2 by the usual procedure of fitting water retention functions (Fig. 2). Thus, parameters w 1 , w 2 , θ r2 , α 2 , and n 2 were estimated using the RETC software (van Genuchten et al., 1991). In addition to θ and h data measured in the dry range in the laboratory (h < − 5000 cm) with a dew point device (Gubiani et al., 2021b), the wet range θ and h data (h ≥ − 200 cm) estimated using equation (2) were used to avoid estimates of w 1 and w 2 that would change the final wet branch provided by equation (2). The limit h = − 200 cm is slightly smaller than the most negative matric potentials measured in the field at the end of the drainage experiments (between − 140 and − 90 cm).
Finally, the inverse solution routine of Hydrus-1D was run once more to find the saturated hydraulic conductivity value (K s2 ) of equation (10) without any additional alteration in all other parameters. Further optimization was also avoided due to the impossibility of evaluating the uncertainty associated with all optimized parameters. This problem is addressed in the discussion section.

Evaluation of results
As the optimization performed by Hydrus-1D is the focus, we numbered the sequence of optimization steps from zero to 5, with optimized parameters shown in parenthesis: step zero (θ s1 , θ r1 , α 1 , n 1 , and K s1 ), step 1 (α 1 and K s1 ), step 2 (θ s1 and θ r1 ), step 3 (α 1 , n 1 , and K s1 ), step 4 (θ s1 and θ r1 again), and step 5 (K s2 ). A trial-and-error procedure was used only in step zero. Note that in steps zero to 4, parameters of the unimodal curves (equations (2) and (3)) were optimized, while for the bimodal curves (equations (6) and (7)), only K s2 was optimized. The standard errors of parameters w 1 , w 2 , θ r2 , α 2 , and n 2 of the bimodal curves fitted with the RETC are presented. The performance of inverse modeling was evaluated using the Root Mean Square Weighted Error (RMSWE) between measured and estimated variables (θ and h), and Akaike (AIC) and Bayesian (BIC) Information Criteria, which asses the adequacy of the model given the increasing number of optimized parameters. Water retention curves for each location and layer were summarized in the "averaged" water retention curve, but all estimated curves are graphically presented.

Results
The sequential optimization gradually reduced the objective function (and RMSWE of measured versus estimated θ and h) to values lower than 0.05 of the initial value at the final step (Fig. 3). Supplementary data (Table S.1) provides all values of RMSWE, AIC, and BIC, for all sequential optimizations carried out in this manuscript.
how sequential optimizations for each location and each drainage experiment gradually reduced the RMSWE value, as well as both information criteria, indicating that each sequential optimization provided a better description of the experimental data as well as a better fitting model. As a result, the estimated θ after the fifth optimization step matched quite well the measured θ during the internal drainage experiments on all dates, locations, and depths (Fig. 4). The greatest deviations between estimated and measured θ were associated with water content variations during drainage. The greatest deviation can be seen for Location 2 (Fig. 4), where the lowest and highest amplitude of θ, considering all layers and dates, were 0.09 and 0.14 cm 3 cm − 3 . A decrease in deviations between estimated and measured θ can be seen at Locations 1, 3, and 4 (Fig. 4), where a wider amplitude of θ occurred. The lowest and highest amplitude of θ was 0.04 and 0.14, 0.17 and 0.22, Fig. 4. Measured (solid lines) and simulated (dotted lines) water contents with calibrated parameters at different depths of four locations (rows) during four drainage experiments (columns). and 0.13 and 0.16 cm 3 cm − 3 , respectively, at Locations 1, 3, and 4. Estimated pressure head also matched quite well their corresponding measurements at the end of drainage (Fig. 5). These data were presented in a scatter plot because only the pressure head measured at the and of drainage were used (see the last paragraph of Section 2.4.1). The highest absolute difference between measured and estimated pressure head was 3 cm.
The estimated parameters (Fig. 6) are located in a physically consistent range, except for some extremely high values of K s2 (6.0E+4, 1.8E+5, 5.4E+5, and 1.4E+6 cm h − 1 ), while all other K s2 values were lower than 3.6E+3 cm h − 1 . The high and apparently inconsistent values of the optimized residual water contents for the unimodal curve (θ r1 had ranged from 0.26 to 0.44 cm 3 cm − 3 ) were reduced to expected values, not higher than 0.12 cm 3 cm − 3 (θ r2 in Fig. 6), when the bimodal function (equation (6)) was used. These values of θ r2 , together with w 2 , α 2 , and n 2 (fitted using RETC, Fig. 6), yielded a final descending segment of the water retention curve (Fig. 7), which is more realistic than the horizontal segment of the unimodal function (Fig. 2) for h lower than − 200 cm (θ between 0.26 and 0.44 cm 3 cm − 3 ). Supplementary data (Table S.2) provide all estimated parameters for each location and each drainage experiment and statistical measures such as Standard Error, tvalue, and 95 % Confidence Limits. The four water retention curves of each layer differed more toward the wet end (Fig. 7), because parameters θ s1 , α 1 , and n 1 of each curve were optimized with independent data, that is, θ and h measured during drainage experiments on November 30, 2020, December 3, 2020, April 25, 2021, and May 22, 2021. On the other hand, the data used in the dry range (θ and h data measured for h < -5000 cm with a dew point device) were the same for all curves. The segment indicating a transition between two different pore systems, which was needed to properly connect the wet and dry segments, is noticeable in all curves (Fig. 7).

Discussion
One of the main problems addressed in this manuscript was to find suitable hydraulic functions for soils with RF. By using θ and h measured during field drainage experiments, were found parameters of unimodal functions (equations (2) and (3)) that yielded simulated water contents quite close to their respective measurements during the entire drainage experiment. However, as discussed in the methodology section (Section 2.4.3), values of θ r1 (between 0.26 and 0.44 cm 3 cm − 3 ) were inconsistent for the coarse-textured soils under investigation (Fig. 1) and would lead to inaccurate simulations of fluxes in drier soils (e.g., root water uptake, infiltration, etc.). We then used the laboratory-measured retention data for h < -5000 cm to correct the dry end of the retention curve and replaced the unimodal functions (equations (2) and (3)) with the bimodal functions (equations (6) and (7)) to keep estimated water contents close to those measured in the wet range.
We can therefore conclude that bimodal functions are more suitable for soils with RF than unimodal functions taking into account the following evidence: (i) bimodal functions yielded simulated water contents quite close to measured ones during the drainage experiments ( Fig. 4), (ii) deviations between simulated and measured θ were lower at the locations with higher contents of RF (1, 3, and 4) than at the location with a smaller content of RF (location 1) (Fig. 4), (iii) optimized parameters were physically consistent except for a few high values of K s2 (Fig. 6), (iv) dry ends of retention curves (Fig. 7) matched well data measured for h < − 5000 cm (from 0.02 to 0.28 cm 3 cm − 3 ), and (v) multimodal functions are commonly used for soils with RF due to the presence of larges pores and cracks surrounding RF, which contrast with textural pores in the aggregated fine material among them (Durner, 1994;Ma and Shao, 2008;Ma et al., 2010;Zhang et al., 2016). As θ and h were not measured for the transition interval between the two porous systems (around h − 200 cm, Fig. 7), it was impossible to evaluate how well the estimated parameters (Fig. 6) describe the real porous system. Future studies can address this issue by adding to the objective function retention data (θ and h) measured beyond the drainage experiment, allowing for water losses due to evaporation and root water uptake when the area is cultivated.
Another challenge addressed in this study was to perform parameter optimization while avoiding ill-posedness and instability problems of the inverse solution, which needs more attention when local gradientbased optimization algorithms, such as the Levenberg-Marquardt scheme, are used (e.g., Š imůnek et al., 2012). Usually, less than five parameters are optimized simultaneously to avoid ill-posedness and instability problems of the inverse solution (Ritter et al., 2003;Abbasi et al., 2003). Our strategy of sequentially optimizing the vector of 15 parameters (θ s1 , θ r1 , α 1 , n 1 , and K s1 ) for three soil layers successfully avoided failure of the Hydrus-1D inverse solution due to non-   6. The range of all optimized parameters at all locations, depths, and times. Only w2 is presented because w1 is complementary to w2 (w1 + w2 = 1, equation (7)).
convergence and gradually reached the minimum of the objective function (RMSWE of measured versus estimated θ and h).
After the initial trial and error approach (step zero), the sequential optimization of α 1 and K s1 (step 1), θ s1 and θ r1 (step 2), α 1 , n 1 , and K s1 (step 3), and again θ s1 and θ r1 (step 4) reduced the RMSWE to < 0.05 of the initial value (Fig. 3) after step 4, which was the last optimization step for parameters of unimodal functions (equations (2) and (3)). In the last step (step 5), in which only K s2 for bimodal functions was optimized using the Hydrus-1D inverse solution, the RMSWE was reduced in several cases even further (Fig. 3). Thus, by combining sequential optimization with bimodal hydraulic functions, in which only a few parameters are optimized while the others are kept constant (Ritter et al., 2003;Abbasi et al., 2003), we successfully performed inverse parameter optimization for layered soils with RF. The lack of data for the transition interval between the two porous systems (from around h = − 200 cm to h = − 5000 cm, Fig. 7), as discussed earlier, prevented the optimization of w 1 , w 2 , θ r2 , α 2 , and n 2 directly using the Hydrus-1D inverse solution and thus these parameters were optimized using RETC. However, even if measured data for this interval were available, fitting all parameters of bimodal functions Fig. 7. Bimodal water retention curves. Gray lines in each graph represent individual retention curves for particular drainage experiments (four dates). Read lines represent the "averaged" curve. Open circles are the measured data with a dew point device.
described by equations (9) and (10) for multilayered soil profiles would remain challenging. A three-layered soil profile would result in a vector with 24 soil hydraulic parameters. It would be unavoidable to fix some parameters to measured values or estimate them independently to deduce the ill-posedness and instability of the inverse algorithm.
Of all parameters, θ s could be measured with relatively low uncertainty. However, θ s cannot be directly measured in soils with RF that prevent the sampling of undisturbed cores. θ s could be estimated as a function of the dry bulk density and the mineral density of the solid fraction when these properties are available. Alternatively, the easiest way of estimating θ s is to assume that its value equals the highest value of θ measured in the soil profile after saturation, which implies ignoring any entrapped air during saturation. An effective Ks could be measured with field methods such as double ring and disc infiltrometers (Baetens et al., 2009;Verbist et al., 2010). However, the effective Ks measured by the infiltration experiment disregards the internal variability of Ks in heterogeneous soil profiles. Furthermore, the Ks as a model input often requires model calibration due to its high sensitivity even when measured values are available (Mertens et al., 2005;Rocha et al., 2006), which implies little usefulness of measured values. Thus, soils with RF and high heterogeneity will probably demand indirect determination of all soil hydraulic parameters.
Differently from earlier studies that evaluated hydraulic proprieties in disturbed samples created with a mix of RF and finer particles collected from structured soils (Ma and Shao, 2008;Beibei et al., 2009;Ma et al., 2010;Parajuli et al., 2017;Naseri et al., 2019), we successfully optimized parameters of hydraulic functions by using inverse modeling of field drainage experiments, in which the real porous system is preserved. Using θ and h data collected in the field (during internal drainage experiments) and laboratory (using a dew point potentiometer), we determined that the combination of the Hydrus-1D inverse solution in sequential mode with parameter estimation using RETC is a suitable strategy to determine hydraulic parameters for complex porous media like soils with RF. These results highlight the suitability of inverse modeling to indirectly determine soil hydraulic properties using transient flow experiments (Šimůnek et al., 1998;Lazarovitch et al., 2007;Jacques et al., 2012;Rashid et al., 2015;Pinheiro et al., 2019;Brunetti et al., 2019). Furthermore, hydraulic parameter estimation by inverse modeling in situ can also contribute to a better understanding of the effects of RF on soil hydrological properties and processes (Zhang et al., 2016).

Conclusions
This study aimed to determine suitable hydraulic functions for soils with rock fragments (RF > 2 mm) and find a sequential optimization strategy to estimate their parameters by inverse modeling using field drainage experiments. Four layered soil profiles with different amounts of RF in different layers (from zero to 99 %) were evaluated in the field using internal drainage experiments and dewpoint potentiometers in the lab. The two main conclusions are: 1-Bimodal functions should be considered to represent hydraulic properties (water retention and conductivity) of soils with RF, i.e., in soils with large pores and cracks surrounding RF, which contrast with textural pores in the aggregated fine material and fine pores inside porous RF. A disadvantage of these functions is an increase in the number of calibrated parameters, leading to ill-posed problems and instability of the optimization algorithm. 2-It is possible to sequentially optimize the van Genuchten parameters of layered soil profiles and successfully avoid failure of the Hydrus-1D inverse solution due to non-convergence. In this study, the parameters of the unimodal van Genuchten-Mualem function (subscript 1) were optimized using the Hydrus-1D inverse solution repeatedly (steps). α 1 and K s1 were optimized in step 1, θ s1 and θ r1 in step 2, α 1 , n 1 , and K s1 in step 3, and θ s1 and θ r1 again in step 4. At the last step, the RMSWE (measured versus estimated water contents and pressure heads) was reduced to less than < 0.02 in most cases. Using the optimized parameters of the unimodal van Genuchten-Mualem function and lab data for pressure heads lower than − 200 cm, all other parameters of the bimodal van Genuchten function (w 1 , w 2 , α 2 , and n 2 ) were estimated by a curve fitting procedure, in which the high incoherent values for θ r1 were also corrected. The Hydrus-1D inverse solution was then run again to optimize only K s2 for the van Genuchten-Mualem bimodal function with all previously estimated parameters constant. In this way, 27 parameters were successfully estimated for layered soils with RF.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.