Estimating macropore parameters for HYDRUS using a meta-model

The vertical change in the number of macropores causes a variation of the relative macroporosity ( w f ) and the effective aggregate width ( d ag ) over the soil profile. Both parameters are used in HYDRUS to represent this variation, increasing the number of parameters and making automated calibration challenging. The working hypothesis is that we can improve an analytical estimation of w f and d ag developed in previous research by inverse estimation with a meta-model for HYDRUS 2D/3D, using disk infiltrometer data of infiltration at zero pressure head. We generate a meta-model that describes the vertical heterogeneity of the macropore number with a general function using four parameters: the relative macroporosity at the soil surface ( w fs ), the effective macropore radius ( r m ), the maximum depth of macropores ( z max ) and the shape parameter of the w f curve ( m ). The meta-model computes the variation of w f and d ag over depth, thus reducing the parameters for automated calibration with HYDRUS. We theoretically described how to directly obtain the meta-model parameters with disk infiltrometer data, providing an example for field conditions. A complete parametrization of matrix and macropore parameters for HYDRUS 2D/3D was generated from these


| INTRODUCTION
The heterogeneous distribution of the macroporosity and the geometry of the macropore network are characterized by unevenly distributed macropores over the soil profile, which results in the number of macropores changing over depth. In Urbina et al. (2019), this situation is mentioned as "heterogeneous macropore geometries." The uneven distribution of macropore numbers over depth was observed in Schon, Mackay, Gray, van Koten, and Dodd (2017) and is probably a common occurrence under field conditions. Simulation of heterogeneous macropore geometries under controlled conditions was performed by Urbina et al. (2019) for HYDRUS-1D (Šimů nek, van Genuchten, & Šejna, 2016) and SWAP (Kroes et al., 2017) computer codes. They found that the variation over depth of the relative macroporosity, w f , and the effective aggregate width, d ag , are critical parameters for simulating this uneven distribution in both models. Hence, advances in the estimation of w f and d ag for field conditions are essential for increasing the accuracy of dual-permeability model simulations.
A method of estimating w f and d ag under field conditions from disk infiltrometer data was developed by Urbina, van Dam, van den Berg, Ritsema, and Tang (2020). The resulting methodology allows an initial estimation of w f and d ag for dual-permeability models such as HYDRUS, SWAP and MACRO (Jarvis & Larsbo, 2012), including different macropore-matrix geometries. The methodology is based on quasi-steadystate one-dimensional water flow conditions. However, steady-state conditions are challenging to achieve for field environments (Šimů nek, Wendroth, & van Genuchten, 1999), and disk infiltrometer infiltration is commonly 3D instead of one-dimensional (Stewart et al., 2016;Alberti & Cey, 2011). Therefore, improvements in Urbina et al.'s (2020) methodology are necessary to adjust w f and d ag estimation for field conditions. The improvement can be performed by calibration using disk infiltrometer data as observations and a 3D dualpermeability model for simulations. This improvement is the primary goal of this research.
Measurements with a disk infiltrometer at zero pressure head allow accounting for infiltration into the macropore and matrix domains simultaneously. The infiltration in both domains can be simulated by HYDRUS 2D/3D, which incorporates axisymmetric flow around the vertical axis and a dual-permeability transfer model (Gerke & van Genuchten, 1993). The simulated data from HYDRUS 2D/3D can be used to update w f , d ag and other macropore parameters by calibration. Prior research applied manual calibration to obtain macropore parameters for HYDRUS 2D/3D using Guelph parameters and disk infiltrometers (Alberti & Cey, 2011;Kodesova, Simunek, Nikodem, & Jirku, 2010). However, fully automated calibration is a better alternative than manual calibration to accurately obtain fitted macropore parameters and uncertainty bounds.
The automated calibration of HYDRUS 2D/3D is challenging because the dual-permeability conceptualization in the model produces more numerical instabilities than the uniform flow conceptualization. Therefore, the ability to find reasonable initial estimates for macropore parameters is critical before calibration, along with reducing model complexity (Arora, Mohanty, & McGuire, 2011). Additionally, the parameters related to water and solute exchange between the macropore and matrix (Gerke & van Genuchten, 1996) vary concerning the density (number) of macropores (Arora et al., 2011;Urbina et al., 2019). Therefore, the number of macropore parameters (w f and d ag ) increases considerably when the number of macropores varies over depth. This situation produces a challenging scenario for the automated calibration of HYDRUS 2D/3D under field conditions because it should be necessary to determine those parameters (w f and d ag ) for each soil layer in the model. A meta-model for HYDRUS 2D/3D can overcome those drawbacks. Metamodels or surrogate models are defined in Asher, Croke, Jakeman, and Peeters (2015) as computationally cheaper models designed to approximate the dominant features of a complex model, which facilitates calibration and uncertainty analysis. An example of meta-modelling is found in Hack-ten Broeke et al. (2016), using the hydrological models SWAP-WOFOST. Ideally, the meta-model parameters should be easy to measure for field conditions and smaller in number than the original model.
We develop a meta-model in the R programming language for reducing the number of macropore parameters (w f and d ag ) in heterogeneous macropore geometries for HYDRUS 2D/3D. The meta-model solves the drawback of the increasing number of parameters in vertical heterogeneity of the macroporosity by a general function. An initial estimation of meta-model parameters with a disk infiltrometer is incorporated, which is next improved by calibration. The working hypothesis is that we can improve the analytical estimation of the relative macroporosity (w f ) and the effective aggregate width (d ag ) developed in Urbina et al. (2020) by inverse estimation with a meta-model for HYDRUS 2D/3D, using disk infiltrometer data of infiltration at zero pressure head. The objectives of this study are as follows: (a) conceptualize a meta-model to reduce the number of dual-permeability parameters during inverse estimation; (b) propose a methodology for an independent measurement of meta-model parameters with disk infiltrometer data; and (c) provide an example of applying the metamodel under field conditions with disk infiltrometer data using HYDRUS 2D/3D for simulations. We refer to HYDRUS-1D/2D/3D collectively as HYDRUS, except when explicitly referring to the dimensionality of the model (e.g., HYDRUS-1D) is necessary.

| HYDRUS model
The HYDRUS code uses the dual-permeability model presented in Gerke and van Genuchten (1993). This model solves two Richards equations, one for the matrix domain and the other for the macropore (fracture) domain, and couples them by a water exchange term: where C is the specific water capacity ( ∂θ ∂h Þ cm −1 , h is the pressure head (cm), K is hydraulic conductivity (cm min −1 ), z is depth (cm) and t is time (min). The subscripts f and m refer to fracture or matrix, respectively. Γ w is the water transfer term (min −1 ) computed as: where β is a dimensionless geometry-dependent coefficient, d ag is the effective aggregate width (cm), which represents an effective distance between the wall of the macropores and the centre of the matrix, γ w is a dimensionless scaling coefficient (usually taken as 0.4) and K a (h) is the effective hydraulic conductivity evaluated in terms of both h f and h m . β can be obtained from d ag and the macropore radius (r m, cm) for cylindrical macropores surrounded by a cylindrical matrix mantle (Gerke & van Genuchten, 1996): The dual-permeability model of Gerke and van Genuchten (1993) was incorporated in HYDRUS 2D/3D, keeping the same conceptualization and parameters mentioned above for the 1D case.
We defined "macropore layers" according to the variation of the macropore numbers in the soil profile, which produce variation in the lateral transfer parameters (w f , d ag, and β) in HYDRUS. Therefore, "macropore layers" is a concept concerning HYDRUS and not related to pedology. Those macropore layers may be in one homogeneous soil matrix horizon described for field conditions. For example, in Figure 1 (III), eleven macropore layers and three matrix layers (or soil horizons) are present. The macropore layers are generated at each macropore ending.
Twenty-two values of w f and d ag are required in HYDRUS to represent the problem in Figure 1 (w f,j, and d ag,j , where j = 1,2,…,11). β is not counted because it is mathematically related to d ag (Equation (3)). The metamodel replaces these twenty-two values with four parameters (next section), which reduces the complexity of modelling heterogeneous macropore geometries. Hence, the meta-model enables HYDRUS to be independent of the number of macropore layers or independent of the effect that produces the vertical variation of the macropore number on lateral transfer parameters (w f , d ag, and β).
One strong assumption in the dual permeability conceptualization of HYDRUS is that the macropore connects the top and bottom boundaries. Therefore, it is not possible to include explicitly dead-end pores. However, dead-end pores are related to an increase in the lateral exchange of water. Thus, they can be implicitly represented by the computation of d ag regarding the vertical distribution of the macropore number. This option was used in Urbina et al. (2019) to represent implicitly deadend pores in controlled laboratory conditions with HYDRUS-1D.

| Meta-model
The meta-model for HYDRUS is based on a continuous decrease of the relative macroporosity (w f ) over the soil profile: where w f_z (cm 3 cm −3 ) is the relative macroporosity at depth z (cm), w fs (cm 3 cm −3 ) is the relative macroporosity at the soil surface (maximum w f ), z is the vertical coordinate (cm; positive upward, zero at the soil surface), z max (cm) is the bottom of the active macropore layer (where w f = 0) and m is a shape factor for the w f curve. The relative macroporosity (w f ) and the effective aggregate width (d ag ) are mathematically related (Urbina et al., 2020). For cylindrical macropores surrounded by a cylindrical matrix mantle d ag is computed as follows (Urbina et al., 2019): The reduction of macropore parameters in HYDRUS with the meta-model is performed by computing w f over depth by Equation (4) and then replacing those values in Equation (5) to obtain d ag over depth. Researchers should note that w f and d ag can be obtained from Equations (4) and (5) with a resolution as small as 1 cm; therefore, several macropore layers can be generated and set in the HYDRUS model. On the other hand, those parameters can be averaged to represent each soil matrix layer (or soil horizon) in HYDRUS (Figure 1 (IV)). The average for macropore parameters introduced in Figure 1 (IV) is generally optional, which means researchers can set in HYDRUS an arbitrary number of macropore layers. The only parameter that must be averaged over the soil profile for HYDRUS is the relative macroporosity (w f ), as is observed in the next sections.
The meta-model is written in R script. The previous equations for the meta-model use input parameters different to those in HYDRUS; therefore, a pre-processing is included in the R script to obtain HYDRUS parameters. This pre-processing included the steps depicted in Figure 1 for averaging the parameters. After the preprocessing, the R script runs HYDRUS internally. The meta-model parameters can be independently estimated under field conditions because they have either physical (w fs , z max, and r m ) or geometrical (m) meaning. In the next section, we propose disk infiltrometers to estimate the meta-model parameters for field conditions.

| Disk infiltrometers
The disk infiltrometer is a device that allows for the measurement of infiltration rates at different pressure heads F I G U R E 1 The schematization for initial measurements of HYDRUS 2D/3D macropore parameters by a meta-model. Part I depicts the variables obtained from disk infiltrometer measurements at two depths, z 1 and z 2 . Part II shows the generation of meta-model parameters from the variables obtained in Part I. Part III illustrates the use of the meta-model under heterogeneous macropore geometries for generating the relative macroporosity w fj , the effective aggregate width d agj and the shape parameter β j for each macropore layer j. Part IV demonstrates how the different parameters are averaged (Avg) in HYDRUS 2D/3D [Color figure can be viewed at wileyonlinelibrary.com] imposed at the base of a disk (Perroux & White, 1988). These pressure heads can be transformed into a macropore radius (Equation (6)) using the Young-Laplace capillarity theory. For water at 20 C and a contact angle of zero between the liquid-vapour-solid interface, the equivalent macropore radius (r m ) for a given pressure head (h) is: where h (cm) is the pressure head imposed at the base of the disk infiltrometer. Let us assume throughout the next sections that a disk infiltrometer is employed at {0, −3, −6 and −10} cm pressure head. We can generate pressure head ranges such as (0, −3), (−3, −6) and (−6, −10) cm, where the macropore water flux is measured at steady-state flow conditions. The minimum radius of each pressure head range is 0.05, 0.025 and 0.015 cm (Equation (6)), respectively. We considered a value of −10 cm for the threshold pressure head between the macropore and matrix domains (Jarvis et al., 2007).
The steady-state flux rates (q M , cm min −1 ) are computed as follows: where q M_0 is the steady-state flux rate measured at zero pressure head and q M[0, −3] is the steady-state flux rate obtained by subtracting q M_3 from q M_0 . A similar computation is performed for q M[−3, −6] and q M[−6, −10] . Dividing Equation (7) by q M_0 yields the flux proportion associated with each macropore size range.
The number of cylindrical macropores per unit area, N (m −2 ), is obtained by applying the Hagen-Poiseuille equation and Darcy's law, assuming laminar flow conditions (Dunn & Phillips, 1991;Watson & Luxmoore, 1986): where η (Kg m −1 s −1 ) is the dynamic viscosity of water, ρ (kg m −3 ) is the density of water and g (m s −2 ) is the gravitational constant. Equation (8) is computed for different pressure heads because of q M and r m (Equation (6)). The total number of macropores, N T (m −2 ), is obtained by: where N [0, −3] is the number of macropores per unit area using q M [0,−3] and r m = 0.05 cm in Equation (5). A similar computation is performed for N [−3, −6] and N [−6, −10] . The relative macroporosity, w f , can be computed using disk infiltrometer data as follows: where A m is the average macropore cross-sectional area (m 2 ). Both N and A m in Equation (10) are computed for different pressure heads. Hence, the final relative macroporosity is obtained as follows: where w f(0, −3) is computed using N [0, −3] and r m = 0.05 cm in Equation (10)

| Computation of meta-model parameters with disk infiltrometers
In previous sections, we indicated how to compute the relative macroporosity with disk infiltrometers, including different pressure head ranges. These concepts are relevant for obtaining initial estimations of meta-model parameters (Equations (4) and (5)). Equation (4) computes the relative macroporosity at depth z in the soil profile. Therefore, to obtain Equation (4) parameters, we need to place the disk infiltrometer in at least two depths, z 1 and z 2 , for different pressure heads (Figure 1 (I)). Researchers should note that by placing the disk infiltrometer at depths z 1 and z 2, an average relative macroporosity is obtained. For example, when placing the disk infiltrometer at the soil surface (z 1 = 0), the measurement of the relative macroporosity by Equation (11) represents an average from z 1 until a depth z (unknown) ( Figure 1_I). The average relative macroporosity and its depth are denoted by the letter "g" in Equations (12) to (15). From previous infiltration rate measurements with the disk infiltrometer at two depths, the relative macroporosity at the soil surface, w fs, and the shape factor parameter, m (see Equation (4)), are obtained as follows: Combining Equations (12) and (13) obtains: where w f g zx represents an average relative macroporosity located at depth g zx (cm) (Equations (12) and (13)). w f g zx is computed by Equation (11); therefore, Equations (14) to Equation (15) can be applied individually to different pressure head ranges.
Equations (14) and (15) allow the estimation of w fs and m in Equation (4). The depths g z1 , g z2 and z max (Equations (12) to (15)) are initially subjectively chosen and are set constant over the different pressure head measurements. However, we anticipate that the depths g z1 and g z2 are updated during calibration.
The parameter d ag can now be computed for different depths and pressure head ranges using Equation (5). However, d ag cannot be set for different pressure head ranges in the HYDRUS model. Therefore, the computation of an average effective aggregate width, d ag (cm), using different pressure head ranges is necessary for obtaining a singular value. This d ag is computed at every macropore layer over depth (Figure 1 (III)). The d ag is obtained from the disk infiltrometer measurements as follows: d ag(0, −3) is computed with Equation (5) where r m = 0.05 cm and w f = w f(0, −3) . The same holds for d ag (−3,−6) and d ag (−6,−10) cm. Equation (16) shows that the average d ag is weighted by water flow.

| Saturated hydraulic conductivity of macropores
The meta-model is not related to the computation of the saturated hydraulic conductivity of the macropores (K sf , cm min −1 ), although K sf is required for the dualpermeability model of HYDRUS. The initial estimation of K sf based on cylindrical macropores can be performed by combining the Hagen-Poiseuille equation and Darcy's law for laminar flow conditions (Dunn & Phillips, 1991;Watson & Luxmoore, 1986): where r m_3 (cm) is the maximum macropore radius computed for pressure head −3 cm (Equation (6)), with similar consideration for r m_6 and r m_10 . The total number of macropores, N T, is computed with Equation (9). Equation (17) for K sf is based on the arithmetic mean of parallel stratified media (Jang, Narsilio, & Santamarina, 2011), weighted by the number of macropores for each pressure head range (e.g., N (0, −3) to N T ). The K sf solution in Equation (17) is valid for horizontal, vertical and sloped macropore directions.
The previous methodology allows us to directly estimate w fs , m, d ag and K sf by disk infiltrometer measurements, whereas g z1 , g z2 and z max are initially subjectively chosen. One way to estimate all these parameters directly without subjective suppositions is by placing the disk infiltrometer within at least five depths (five unknowns in Equations (12) to (15)). However, that procedure is out of the scope of this research.

| Field example
A field site with Spruce trees was selected to test the meta-model developed in previous sections. The field site is in Valkenburg, the Netherlands (N50 51.515/ E005 53.684, elevation 117 m). The 20-year average precipitation and average temperature until 2019 are 763.3 mm and 10.84 C, whereas the 10-year average precipitation and average temperature until 2019 are 726.4 mm and 10.85 C. The field site was not ploughed for 8 years. Therefore, a significant macropore presence was expected. Only cylindrical macropores were observed in the field. The textural class changed from clay loam at the topsoil to clay at the subsoil. The soil was found to be hydrophilic at the top 50 cm; below that, no measurements were performed.
Undisturbed soil core samples of 8 cm high and 10 cm in diameter were taken at the soil surface, and at 18.5 and 33-cm depth, after the disk infiltrometer measurement. The measurements were obtained directly below the place where the disk infiltrometer was set. The samples' depths were selected to represent each soil horizon present to a depth of 50 cm.
The undisturbed soil cores were analysed in the Soil Hydro Physics Laboratory of Wageningen University & Research Centre. The soil hydraulic matrix parameters were determined using the evaporation method. Four mini tensiometers were installed at 2-cm intervals, and at 1 cm from the top and bottom boundaries of the soil core.
The measured soil water pressure head and soil evaporation were analysed through calibration using HYDRUS-1D (Šimů nek, van Genuchten, & Wendroth, 1998), setting the residual water content (θ r , cm 3 cm −3 ) to the default value provided by HYDRUS-1D for clay. Next, the saturated water content (θ s , cm 3 cm −3 ), the inverse of the air entry value (α, cm −1 ), the pore size distribution (n), the saturated hydraulic conductivity (K s , cm min −1 ) and the pore connectivity (l) parameters of the matrix were calibrated. Care was taken to ensure the same measured and simulated evaporation at the end of the optimization.
The derived van Genuchten-Mualem parameters (θ r , θ s , α, n, K s , l) and the bulk density (D a , g cm −3 ) are presented in Table 1. The bulk density was obtained using the undisturbed soil cores at the end of the evaporation method (ISO 11272, 1993).
A tension disk infiltrometer of 20 cm diameter (Commonwealth Scientific and Industrial Research Organisation, Disk Permeameter) was set at z 1 = 0 and z 2 = 30 cm depth. For each measurement, around 1 cm of coarse sand was placed and levelled. At the soil surface, all grass was removed by cutting it as much as possible.
Infiltration rate measurements were performed at the soil surface (z 1 ) for 0, −3 and − 6 cm pressure head until a steady-state was reached. The measurement at 0 cm pressure head was performed first. One hour later, measurements at −6 and −3 cm pressure head were performed. After that, a pit was gently dug directly below the position where the previous measurement was taken. The disk infiltrometer was placed in this new position (z 2 ) , and the infiltration rate was measured only for 0 cm pressure head at steady-state conditions.
The infiltration rate measured at the soil surface for 0 cm pressure head was measured in short time intervals until steady-state because that information is used for calibration with HYDRUS 2D/3D. A total of 34 measurements were performed in 58 min with a disk infiltrometer.
The steady-state flow rates, q M_0 , q M [0,−3] and q M [−3,−6] in Equation (7), were directly obtained from measurements at depth z 1 = 0 cm. The steady-state flow rate, q M [−6,−10], was estimated for depth z 1 as follows: The steady-state flow rate for the matrix, q M[matrix] , was estimated as follows (Moene & van Dam, 2013): where N h = 3 is the number of matrix layers, L mk (cm) is the thickness of each matrix layer k, and K sk (cm min −1 ) is the saturated hydraulic conductivity of the matrix layer k. We assumed that q M[matrix] is equal to the composite saturated hydraulic conductivity obtained by Equation (19) under steady-state conditions (Watson & Luxmoore, 1986). The number of matrix layers, L mk, and the saturated hydraulic conductivity, K sk of each matrix layer (k), are found in Table 1. The measurement at depth z 2 = 30 cm was only performed for a steady-state flow rate at 0 cm pressure head, q M_0. The ratio between q M_0 measured at depth z 2 and z 1 was 0.39 ( Table 2). The ratio was multiplied with q M_3 , q M_6 and q M_10 (Equation (7)) measured and estimated at the soil surface to generate those values for z 2 . Therefore, we assume the decrease of steady-state flow, q M_0 in z 2 relative to z 1, proportionally reduced q M_3 , q M_6 and q M_10 . A summary of the steady-state flow rate measurements is given in Table 2.
The meta-model parameters w fs and m (Equation (12) to Equation (15), respectively) were computed by setting the depths g z1 = −15, g z2 = −50 and z max = −100 cm. The values of w f_gz1 = 9.14 * 10 −4 and w f_gz2 = 3.26 * 10 −4 were computed using Equations (12) and (13) (See Supplementary Material). The outcomes were w fs = 1.25 * 10 −3 and m = 1.94. The values of g z1 and g z2 were chosen to be close to the measurement depth of the disk infiltrometer (z 1 and z 2 ) because we expected many macropores close to the soil surface. The value of z max means that we do not expect macropores below 100 cm depth.
The initial gravimetric soil water content was measured at 45 and 40 cm away from the outer disk infiltrometer ring, and both samples were merged into T A B L E 1 The van Genuchten-Mualem parameters (θ r , θ s , α, n, K s , l) and the bulk density (D a ) for the three soil horizons described in the field site one sample. The initial water content was measured at every 10 cm up to a maximum depth of 80 cm. The samples were weighed on-site on a small scale and then oven-dried at 105 C. The initial gravimetric soil water content was transformed into volumetric water content with the bulk density of the soil ( Figure S2, Supplementary Material).

| HYDRUS 2D/3D set up
The axisymmetric flow around the vertical axis was set in HYDRUS 2D/3D to simulate disk infiltrometer infiltration. The disk infiltrometer radius was represented by a red line of 10 cm at the top left corner (Figure 2). A square domain of 50 × 50 cm was set in HYDRUS 2D/3D, and three matrix and macropore layers were included in the model (Figure 2). The matrix and macropore layers corresponded to the soil horizons observed under field conditions. The dual-permeability module of HYDRUS 2D/3D was applied for the simulations. The time and space discretization for HYDRUS 2D/3D were set to obtain a mass balance error lower than 0.5% at every time step. The initial soil water content was set as a pressure head, transforming the initial volumetric water content ( Figure S2, Supplementary Material) using the soil hydraulic matrix parameters ( Table 1). The boundary condition for the disk infiltrometer (red line in Figure 2) was set as a constant head equal to zero. The bottom boundary condition was set as free drainage. All other boundaries were set as no flux boundaries. The model was run for 58 min, which corresponds to the disk infiltrometer measurement at zero pressure head in the soil surface.
The previous estimations of meta-model parameters (w fs , m and z max ), along with Table 2, are sufficient to generate initial estimations of w f, d ag and K sf in HYDRUS 2D/3D. Researchers should notice that the meta-model was generated at centimetre-scale until z max = 100 cm (see Table S1, Supplementary Material). The simulation with HYDRUS 2D/3D terminates at 50-cm depth simply because, in this example, we did not measure the soil matrix at deeper layers.
The relative macroporosity cannot vary over depth in the HYDRUS 2D/3D version 2. xx. Therefore, an average value over depth should be computed before incorporating that parameter into the model (see Urbina et al., 2019). This issue is depicted in Figure 1 (IV), where the 11 w fj are averaged to generate the input parameter for HYDRUS 2D/3D. It was decided to average d ag and β to generate a unique value for each soil matrix layer (Figure 1 (IV)). Two values of K sf were obtained, applying Equation (14) from disk infiltrometer measurements at z 1 and z 2 (Figure 1 (I)). In this research, the average of K sf at the two depths is computed and set constant over depth for HYDRUS 2D/3D (Figure 1 (IV)). However, this is optional because the HYDRUS code allows different K sf over depth. The above-mentioned computations of w fs , m, w f , d ag , β and K sf (Table 3) are available in the Supplementary Material.
The remaining macropore parameters in Table 3 were estimated as follows. The van Genuchten macropore T A B L E 2 The steady-state flux rate measured at 0 (q M_0 ), −3 (q M_3 ) and −6 (q M_6 ) cm pressure head  (18)  parameters (subscript f ) θ rf and θ sf were set to represent a macropore that dries completely and transport water almost using the whole pore volume. The parameters α f and n f were set like the ones for sandy soil. The pore connectivity parameter of macropores l f and the empirical factor γ (−) were fixed to 0.5 and 0.4, respectively. The effective hydraulic conductivity of the fracture-matrix interface, K a (cm min −1 ), was set equal to the saturated hydraulic conductivity of the matrix, K s (Table 1).

| Calibration strategy
The model-independent Parameter Estimation δ Uncertainty Analysis (PEST) package (Doherty, 2015) was used for the calibration of the initially estimated meta-model and HYDRUS 2D/3D parameters ( Table 3). The strategy was to calibrate the meta-model parameters running the R script with PEST. Therefore, PEST reads the input and output files generated by the R script for the meta-model, which runs HYDRUS 2D/3D internally after pre-and post-processing. This calibration was made possible by transforming the R script into a Batch file for Windows. The meta-model parameter selected for calibration was the relative macroporosity at the soil surface, w fs , because it can generate other model parameters via the analytical equations presented in the previous sections. Therefore, the change in w fs performed by PEST at each iteration produces a variation in other model parameters. This interaction between parameters during calibration is illustrated by a tree chart (Figure 3).
The red boxes (Figure 3) mean that those variables were fixed during calibration (e.g., w fs[0,-3] ). The reason for fixed w fs[0,-3] is related to an overestimation of the water flow by the physical assumptions in the methodology (e.g., Hagen-Poiseuille equation, see Urbina et al., 2020). The overestimation of water flow produces an underestimation of w fs and an overestimation of K sf . One way to increase the final value of w fs and decrease the value of K fs simultaneously is by fixing w fs[0,-3] during calibration.
PEST was operated in "estimation mode," which minimizes the objective function with the Gauss-Marquardt-Levenberg method. The observed data (cumulative infiltration at 0 cm pressure head in the soil surface) were compared against the simulated cumulative infiltration of HYDRUS 2D/3D at 0 cm pressure head (red line in Figure 2). The observed data were divided into three groups, which facilitates the calculation of the weights. The PEST module PWTADJ1 computed the weights of each observational data group. They have been calculated with the condition that all observation data groups had the same importance in the objective function (ϕ). Therefore, the three cumulative infiltration groups had the same importance after the residual was multiplied by the weight (Equation (20)): where Nt is the total number of observations, w i is the weight associated with the iʼth observation, and r ei (cm) is the iʼth residual (difference between model output and measurement).
T A B L E 3 Initial estimation of macropore parameters for HYDRUS 2D/3D at each soil horizon

| RESULTS
The relative macroporosity at the soil surface, w fs (Equation (4)), was the only parameter changed in calibration, increasing~3.5 times relative to its initial measurement. This increase in w fs means that the saturated hydraulic conductivity of the macropores, K sf , and the effective aggregate width, d ag , decreased relative to their initial measurements (Figure 3 and Table 4). The decrease in K sf is explained because we only allow changing w fs[−3,-6] and w fs[−6,-10] during calibration; therefore, the number of smaller diameter macropores increases (see Equation (17)). The physical meaning of d ag is an average distance between the macropore wall and the matrix centre. An increase in w fs after calibration means a higher number of macropores that reduce this distance or d ag . The uncertainty bounds were narrow, indicating a good estimation by PEST (Table 4). The other parameters listed in Table 4 were pegged to w fs during calibration. Consequently, those parameters do not have uncertainty bounds (Table 4). The cumulative infiltration simulated by HYDRUS 2D/3D from the initial measurements (Table 3) was lower than that measured under field conditions ( Figure 4). The Nash Sutcliffe coefficient for that simulation was −0.09. The Nash Sutcliffe coefficient for cumulative infiltration simulated by HYDRUS 2D/3D after the calibration of w fs was 0.88.
The increase of w fs after calibration produces an adjustment in the position of the w f average (located at depth, g zx ; see Equations (12) to (15)) ( Figure 5). The subjectively chosen values of the depths g z1 = −15 cm and g z2 = −50 cm change after calibration to g z1 = −56 and g z2 = −74 cm. Therefore, calibration of w fs updates the arbitrarily chosen depths g zx . Recall that the average relative macroporosity at each depth, g z1 and g z2 , is computed from the w f curve generated by the meta-model with Equation (4).

| Initial measurements of metamodel and macropore parameters
Meta-model and macropore parameters (w fs , m, K sf , w f , d ag and β) were directly obtained over the relevant soil depth from disk infiltrometer measurements, whereas other metamodel and macropore parameters (g z1 , g z2 , z max , θ rf , θ sf , α f , n f , γ, l f and K a ) were initially guesses that were later updated by calibration (g z1 , g z2 ) or based on previous research. Köhne, Mohanty, and Šimů nek (2006)  Simunek, Jarvis, and van Genuchten (2006) used similar assumptions for the setting of α f , n f , γ and l f with HYDRUS 2D/3D under field conditions. Köhne and Mohanty (2005) and Gerke and van Genuchten (1993) considered the initial setup of K a equal to the matrix K s as performed in this research. Therefore, the initial parametrization performed for HYDRUS 2D/3D for the above parameters is concordant with previous research. We report that HYDRUS 2D/3D simulations were stable. Therefore, the initial macropore parameters previously obtained by the meta-model and previous research (Table 3) are sound for the dual-permeability concept of HYDRUS 2D/3D. This finding is relevant because researchers can apply similar assumptions for setting parameters in HYDRUS 2D/3D.
The average effective macropore width, d ag , and K sf were weighted proportionally by water flow or number proportion in Equations (16) and (17). We do not know if that is the best way to weight d ag and K sf to obtain an accurate initial measurement. Therefore, this is an open topic for future research.
In this research, we only measured the soil hydraulic parameters of the matrix for the top 50 cm of soil. That information was set for the HYDRUS 2D/3D model with a square domain of 50 × 50 cm and free drainage bottom boundary condition. This implies that the model cannot account for a drastic reduction in infiltration rate when the macropore domain is filled with water.

| Calibration
The Nash Sutcliffe coefficient obtained from the initial measurements of the meta-model and macropore parameters by disk infiltrometer data was −0.09. Therefore, calibration was necessary and improved the Nash Sutcliffe coefficient substantially to 0.88. A visual inspection of Figure 4 indicates a mediocre agreement between the calibrated and observed cumulative infiltration curve at the beginning, where a strong capillarity effect is seen. The effect of capillarity forces could have been produced because the sand material utilized was perhaps not fully saturated or simply because of unsaturated matrix flow at the beginning of infiltration. HYDRUS 2D/3D could not reproduce that initial non-linear curve in Figure 4 because we only calibrate one meta-model parameter (w fs ) related to macropore flow, with a negligible effect on matrix flow. Additionally, the matrix parameters directly estimated under laboratory conditions by the evaporation method are not error free. Matrix parameters can be calibrated under field conditions using cumulative infiltration from disk infiltrometer data at negative suctions (Kodesova et al., 2010).
The calibration of w fs was tied to other macropore model parameters. Our prior information indicated that water flow is overestimated by the physical assumptions, thereby underestimating w fs and overestimating K sf . The increase of w fs~3 .5 times after automated calibration of the initial value confirms our prior information. Alberti and Cey (2011) utilized the Hagen-Poiseuille equation for the initial measurement of K sf. They also concluded that the initial value is overestimated and should be decreased by manual calibration. The overestimation of macropore flow by the Hagen-Poiseuille equation was previously discussed in Dunn and Phillips (1991) and Watson and Luxmoore (1986).
The final relative macroporosity at the soil surface (w fs ) obtained by calibration might still be underestimated. The value of K sf is still high compared to other studies (Arora et al., 2011or Urbina et al., 2019. However, it is comparable to the one obtained in Alberti and Cey (2011). We think that the high initial value of K sf obtained by Equation (17) influences the calibration of w fs . In other words, starting the inverse estimation with a lower value of K sf would increase the final w fs value obtained by calibration more than~3.5 times. We performed a calibration of K sf , keeping constant w fs to its initial value (0.0013 cm 3 cm −3 ; Table 4). The final calibrated value of K sf was 699.2 cm min −1 with 95% uncertainty bounds (622.8-784.9) and Nash Sutcliffe coefficient of 0.88. Therefore, a good match between observation and simulations of the cumulative infiltration curve at zero pressure head can be obtained by either increasing w fs or K sf . This indicates that an accurate initial estimation of K sf is necessary to estimate w fs by calibration. In future research, Equation (17) can be modified to include filmflow (Nimmo, 2010) or rivulet flow (Germann et al., F I G U R E 5 Relative macroporosity over depth computed by the meta-model from initial measurements (wf_ini) and after calibration of w fs (wf_cal). w f_gz1 and w f_gz2 are the average relative macroporosity computed at depths z 1 = 0 and z 2 = 30 cm, respectively, by disk infiltrometer. The subjectively chosen values g z1 = −15 cm and g z2 = −50 cm change after calibration to g z1 = −56 and g z2 = −74 cm [Color figure can be viewed at wileyonlinelibrary.com] 2007) for macropore diameters higher than 3 mm (Germann, 1987). The previous adjustments would reduce the water flux density and thus reduce the value of K sf compared with the one obtained by Equation (17).
Future field applications should also consider the storage change to improve the calibration with disk infiltrometer data (Šimů nek et al., 1999). The storage change is measured by subtracting the final soil water content after disk infiltrometer infiltration from the initial soil water content. Including the storage change over depth and the measurement of cumulative infiltration at 0 cm pressure head at the soil surface introduce additional information for the calibration of w fs and perhaps would allow calibration of an additional parameter.
Previous studies with HYDRUS advised lumping the parameters related to lateral flow, such as the effective hydraulic conductivity of the fracture-matrix interface (K a ), the scaling factor (γ), the macropore shape (β) and the effective aggregate width (d ag ), into one single parameter during calibration (Haws, Rao, Simunek, & Poyer, 2005). Although that solution is practical, it is not easy to provide an initial estimation for the lumped parameter. Additionally, that parameter remains constant over depth; therefore, it cannot represent w f and d ag variation in heterogeneous macropore geometries. Instead, the meta-model presented here allows finding initial measurements of meta-model parameters by disk infiltrometer data. Moreover, during calibration, the parameters can be tied to other macropore parameters because they are mathematically related.
In this research, we observed only cylindrical macropores for field conditions. The methodology developed in Urbina et al. (2020) details the measurement of w f for different macropore-matrix geometries such as rings, hexagons, bricks and rectangular slabs using disk infiltrometer data. Therefore, the initial estimation of meta-model parameters can be used for different macropore-matrix geometries, increasing the accuracy of the initial estimation of w f and d ag .
The feasibility of this research should be analysed in terms of the meta-model and the initial estimation of meta-model parameters by disk infiltrometers. The metamodel included in Equations (4) and (5) and Figure 1 can be used in any inverse estimation for laboratory or field conditions when the number of macropores variates over depth. The initial estimation of meta-model parameters by the disk infiltrometer (Urbina et al., 2020) and its posterior update by inverse estimation (this research) should be used as an initial guess of meta-model parameters for detailed plot studies or risk assessment studies where data for calibration are scarce.

| CONCLUSIONS
A meta-model was developed for HYDRUS 2D/3D in the R programming language. The meta-model reduces the number of macropore dual-permeability model parameters needed to simulate field conditions with a vertical variation of the macropore number. The meta-model computes the variation of the relative macroporosity, w f , and the effective aggregate width, d ag , over depth using only four parameters. The new parameters are the relative macroporosity at the soil surface (w fs ), the effective macropore radius (r m ), the maximum depth of macropores (z max ) and a shape parameter for the w f curve (m). All the meta-model parameters have either physical or geometrical meaning; therefore, they can be independently estimated for field conditions. Independent estimation of meta-model parameters is theorized using disk infiltrometer data. The methodology requires data from disk infiltrometer measurements, at least within two depths with different pressure heads, including the 0 cm pressure head.
The previous methodology was applied for field conditions using a disk infiltrometer at two depths for different pressure heads. The initial parameter values of HYDRUS 2D/3D were improved by calibration with the cumulative infiltration data obtained from the disk infiltrometer at the soil surface for a 0 cm pressure head. Only w fs was calibrated because the physical assumptions of the methodology underestimate this parameter and overestimate K sf . The calibration was performed with PEST and was stable. The calibrated value of w fs was~3.5 times higher than the initial measurement, which is in line with our prior information. We show how the calibration of w fs influenced the values of other macropore parameters, as they are mathematically related. We consider these contributions useful for the initial setup of dual-permeability models to simulate water flow and contaminant leaching in regional risk assessments and detailed plot studies.