Inverse modelling for predicting both water and nitrate movement in a structured-clay soil (Red Ferrosol)

Horizontal solute absorption

Absorption of a NO₃^– solution by the soil was measured in horizontal columns between 17 and 50 cm in length depending on absorption periods. The air dry soil was moistened before packing into the columns. Columns were packed (using a drop hammer) with relatively dry soil (water content 0.15 g g⁻¹) in two to three g increments to achieve a bulk density of close to 1.03 g cm⁻³, which is similar to that measured in the field. Using a Mariotte bottle, a 110 µmol_c cm⁻³ NO₃^– solution was applied to the inlet of the soil column at zero suction. The outlet of the column remained open, tamped with cotton wool to hold the soil in place. Flow was stopped at set times and the column divided into sections that ranged from 1 to 2.5 cm. Short sections were place near the wetting front to provide an accurate measure of the solute and water contents in this area. The soil sections were transferred to tubes and weighed to determine moist weight.

Two types of column experiments were conducted, and are referred to as Set A and Set B. Each individual experiment in both Set A and Set B used a freshly prepared soil column. Set A consisted of five experiments with infiltration times of 26, 30, 43, 47, and 70 min. In these experiments, water-soluble NO₃^– vs distance in the column was determined by adding deionised water to each column section to make a soil-to-water ratio of 1:5.5 (SD ± 0.4). The soil plus water was weighed. Samples were shaken for 4 h in an end-over-end shaker, centrifuged at 9,800 m s⁻² for 10 min, and the supernatant decanted and weighed. The soil remaining in the tubes was also weighed.

Set B involved four columns with absorption times for two of 80 and 320 min for the others. Duplicate columns in this series were either (i) extracted as in Set A or (ii) the soil solution was extracted using the TFE method described by Phillips & Bond (1989). The adsorbed NO₃^– was extracted by adding a volume of 2M KCl to form a 1:5.5 (SD ± 0.5) soil:solution ratio (Rayment & Higginson, 1992).

The tubes were reweighed and shaken for 1 h to extract adsorbed NO₃^–. The soil plus 2M KCl samples were centrifuged at 9,800 m s⁻² for 10 min, and the supernatant decanted in preweighed falcon tubes and weighed. The soil remaining in the tubes was washed twice by shaking for 30 min in 20 cm³ deionised water to remove residual salts. The tubes were centrifuged at 9,800 m s⁻² for 10 min and the wash solution discarded. The washed soil was oven-dried at 105 °C and weighed to give the oven-dry mass of soil in each section. Water and KCl extracts were analysed for NO₃^––N on an Alpkem autoanalyser (Alpkem, 1992).

Point-source solute infiltration

Perspex wedges were constructed with the same dimensions described by Li, Zhang & Ren (2003) and packed with dry soil (water content of 0.2 g g⁻¹) to a bulk density of 0.95 g cm⁻³. All solutions were applied to the 15° corner of the wedge at a depth of five cm (Fig. 1) with a peristaltic pump set to deliver solution at 50 cm³ h⁻¹, equivalent to a dripper output of 1,200 cm³ h⁻¹ in a 360° flow environment.

Figure 2 shows the two irrigation scenarios applied to the wedge experiments. Both treatments were irrigated for 0.5 h (equivalent to 25 cm³ of solution) with the NO₃^– solution applied to the horizontal columns, that is 110 µmol_c NO₃^– cm⁻³. This was immediately followed by a 1.5 h application of solute free water (75 cm³). The soil was either (i) sampled immediately after the water application (Scenario A) or (ii) allowed to rest for 16 h before irrigating again with water for 6 h (300 cm³) prior to sampling (Scenario B). That is, a total of 375 cm³ of water was applied to the wedge in Scenario B with samples being taken 24 h after initial application of the solute. A single replicate was used for Scenario A and Scenario B was duplicated.

Soil was sampled using the method described by Li, Zhang & Ren (2003). Briefly, a five cm grid was placed over the column and a soil core (two cm internal diameter) taken from the centre of each grid. Additional soil samples were taken at the edge of the wetting front. The soil from the core was sub-sampled to determine gravimetric water content and total NO₃^– concentration (solution and adsorbed). Wet soil was extracted with 2M KCl (1:10 soil: KCl ratio) and analysed on an Alpkem autoanalyser (Rayment & Higginson, 1992). Calculations of total NO₃^– partitioned into the solution and the amount of adsorbed NO₃^– were done using the adsorption isotherm determined from the horizontal column data (described below).

Nitrate adsorption isotherm

Partitioning of NO₃^– between the adsorbed and solution phases, over the range of soil solution concentrations in the columns, was measured by displacing the soil solution with TFE (Phillips & Bond, 1989). The Langmuir equation (Eq. (1)) was then fitted to the data to describe the NO₃^– adsorption isotherm.

where C_a is the concentration of adsorbed solute (µmol_c g⁻¹), C_w is the concentration of solute in the soil solution (µmol_c cm⁻³), C_max is the maximum amount of solute that can be adsorbed by the soil (g cm⁻³), and ϕ determines the magnitude of the initial slope of the isotherm (Sposito, 1989).

C_max and ϕ were determined using the Gauss–Newton nonlinear curve model in the statistical program SAS (version 9.1; SAS, Cary, NC, USA). C_max and ϕ were determined to be 23.17 (95% confidence interval (CI) ± 3.43) and 0.00766 (95% CI ± 0.00194), respectively.

Hydrus flow equations

Water flow is described by the Richards equation modified to describe horizontal flow in one dimension with no loss of water due to evaporation of root uptake (Šimůnek et al., 2008): $$\frac{\partial \text{θ}}{\partial t}=\frac{\partial }{\partial x}\left(K\frac{\partial \text{ψ}}{\partial x}\right),$$ where θ is the volumetric water content (cm³ cm⁻³), t is time (min), x is the horizontal distance (cm), ψ is the water tension (cm), and K is the hydraulic conductivity (cm min⁻¹) given by: $$K(\text{ψ},x)=K_{\text{sat}}(x)K_r(h,x),$$ where K_r is the relative hydraulic conductivity (no unit) and K_sat is the saturated hydraulic conductivity (cm min⁻¹; Šimůnek et al., 2008).

The modified form of the Richards equation that describes water movement in two dimensions assuming no loss of water through root uptake or evaporation can be written (Šimůnek, Van Genuchten & Šejna, 2006): $$\frac{\partial \text{θ}}{\partial t}=\frac{\partial }{\partial x_i}\left[K\left(K_{ij}^A\frac{\partial \text{ψ}}{\partial x_j}\right)+K_{iz}^A\right]$$ where x_i (i = 1,2) are the spatial coordinates (cm), and K_ij^A and K_iz^A are components of a dimensionless anisotropy tensor K^A. Assuming flow is isotropic (i.e. K is equal in horizontal and vertical directions) the diagonal entries of K_ij^A equal one and the off-diagonal entries equal zero (Šimůnek, Van Genuchten & Šejna, 2006).

In the two-dimensional flow scenario K is given by: $$K(\text{ψ},x)=K_{\text{sat}}(x,z)K_r(h,x,z).$$

If the modified form of the Richards equation is applied to planar flow in a vertical cross section, x₁ = x is the horizontal coordinate and x₂ = z is the vertical coordinate. This equation can also describe axisymmetric flow when x₁ = x represents a radial coordinate (Gardenas et al., 2005). The transcripts i and j denote either the x or z coordinate.

To solve the Richards equation, Hydrus implements the soil hydraulic functions of the van Genuchten–Mualem to describe unsaturated hydraulic conductivity in terms of soil water retention parameters (Šimůnek, Van Genuchten & Šejna, 2006). Water retention is described by Šimůnek, Van Genuchten & Šejna (2006) as: $$\text{θ}(\text{ψ})=\{\begin{array}{l}{\text{θ}}_r\frac{{\text{θ}}_s-{\text{θ}}_r}{{[1+{\left|\text{αψ}\right|}^n]}^m}\text{\hspace{1em}\hspace{1em}\hspace{1em}ψ}<0\\ {\text{θ}}_s\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}ψ}\ge 0\end{array}$$ and unsaturated conductivity is written (Šimůnek, Van Genuchten & Šejna, 2006): $$K(\text{ψ})=K_s{S}_e^l{\left[1-{(1-S_e^{1/m})}^m\right]}^2\text{\hspace{1em}\hspace{1em}\hspace{1em}ψ}<0$$ where: $$m=1-1/n,n>1$$ and: $$S_e=\frac{\text{θ}-{\text{θ}}_r}{{\text{θ}}_s-{\text{θ}}_r}$$ In the above equations θ_r is the residual water content (cm³ cm⁻³), θ_s is the saturated water content (cm³ cm⁻³), α (cm⁻¹), n (no unit), and l (no unit) are curve fitting parameters for the hydraulic conductivity function and S_e is the effective water content (cm³ cm⁻³).

When the van Genuchten–Mualem model is used to solve Richards equation in Hydrus there are six soil hydraulic parameters required (θ_r, θ_s, α, n, l, and K_sat).

The Langmuir equation was used to predict the reactive solute transport. In Hydrus, desorption of a nontransforming solute is described by the generalised nonlinear equation (Šimůnek, Van Genuchten & Šejna, 2006): (2) $$C_a=\frac{k_s{C}_w^{\text{ω}}}{1+\text{ϕ}{C}_w^{\text{ω}}},$$ and: (3) $$\frac{\partial C_a}{\partial t}=\frac{k_s\text{ω}{C}_w^{\text{ω}-1}}{{(1+\text{ϕ}{C}_w^{\text{ω}})}^2}\frac{\partial C_w}{\partial t}+\frac{C_w^{\text{ω}}}{1+\text{ϕ}{C}_w^{\text{ω}}}\frac{\partial k_s}{\partial t}-\frac{k_s{C}_w^{2\text{ω}}}{{(1+\text{ϕ}{C}_w^{\text{ω}})}^2}\frac{\partial \text{ϕ}}{\partial t}+\frac{k_s{C}_w^{\text{ω}}\ln C_w}{{(1+\text{ϕ}{C}_w^{\text{ω}})}^2}\frac{\partial \text{ω}}{\partial t},$$ where k_s (cm³ g⁻¹), ω (dimensionless), and ϕ (cm³ g⁻¹) are constants. In the case of the Langmuir equation, k_s = C_max ϕ (where C_max and ϕ are the Langmuir equation constants from Eq. (1)) and ω = 1.

Estimation of flow equations parameters

Soil hydraulic parameters for the van Genuchten equation were determined using the inverse optimisation procedure in Hydrus 1D (Šimůnek et al., 2008). Water profile data from the soil columns in Set A were used, after removing obvious outliers in the data that were determined to be due to soil loss during column sampling (Hopmans et al., 2002). The initial water content was set to 0.15 (cm³ cm⁻³), the measured water content of the repacked columns, for all optimisations. Free water absorption was simulated by applying constant water content boundary condition of 0.64 (cm³ cm⁻³) to the opening of the column. The lower boundary condition was set to free drainage. The maximum number of iterations was set to 50 and 86 water content points across five time steps were used in the inverse scenario from the Set A column experiments. Weighting of inverse data was by standard deviation and an equal weighting was applied to all the water content values. Residual soil water content (θ_r), θ_s, and K_sat were set or optimised depending on the optimisation scenario (Fit All is the term applied when all parameters were optimised and Set Measured is used when θ_r, θ_s, and K_sat were set to independently measured values). The remaining empirical parameters, α, n, and l were fitted by running the inverse parameter estimation option in Hydrus 1D. The initial values of α, n, and l were based on the default values given for the Loam soil in Hydrus 1D (α = 1.56, n = 0.173, and l = 0.5).

Initial estimates for θ_r and θ_s were 0.05 and 0.58 (cm³ cm⁻³) and K_sat was set to 0.10 cm min⁻¹. Saturated soil water content (θ_s) and K_sat values were independently measured in falling head K_sat experiments (Reynolds et al., 2002). A column of water (4.2 cm in diameter and 12 cm high) was applied to wet repacked soil core packed to a bulk density of 1.0 with air dry soil sieved to <2 mm. The soil cores were two cm high and had an internal diameter the same as the water column sitting above it. Triplicate measurements of conductivity were recorded on four separate cores. Residual soil water content (θ_r) was estimated based on the air-dry soil water content.

The Rosetta PTF model (Schaap, Leij & Van Genuchten, 2001) was used to estimate soil hydraulic properties as a comparison against the inverse modelling method. The soil particle size measurements (Table 1) and bulk density (1.03 g cm⁻³) were used to estimate soil hydraulic parameters. In a second prediction, moisture retention at −33 and −1,500 kPa (0.34 and 0.22 cm³ cm⁻³, respectively) were also included as inputs into Rosetta. The water content at −33 kPa was determined on a suction table apparatus and −1,500 kPa were determined using pressure plate apparatus (Cresswell, 2002).

Modelling water and solute absorption in soil wedges

Hydrus 2D/3D was used to model water distribution in the horizontal column experiments based on the same initial and boundary conditions used in Hydrus 1D during inverse optimisation. Horizontal flow in Hydrus 2D/3D was simulated by setting the geometry to a 2D horizontal plane. The geometry of the flow domain was set to a column two cm in diameter and 50 cm long. Soil hydraulic parameters having the lowest values for the objective function were used to simulate water absorption. Nitrate absorption was predicted by applying a third-type (Cauchy) solute boundary at the inlet of the column (x = 0 cm) at a constant concentration of 110 µmol NO₃^– cm⁻³. Bulk density was set to the measured value of 1.03 g cm⁻³, longitudinal and transverse dispersivities were set to 0.3 and 0.03 cm, respectively (Ajdary et al., 2007), and the diffusion coefficient was neglected as it was considered negligible relative to the dispersion (Hanson, Šimůnek & Hopmans, 2006; Ajdary et al., 2007). Parameters for the Langmuir equation (Eq. (1)) to describe NO₃^– adsorption were determined from the data measured in column Set B.

Modelling water and solute in horizontal columns

To simulate water and NO₃^– distribution in the wedge experiments, a 40 × 40 cm flow domain was created in a 2D axisymmetrical vertical flow geometry. The infiltration point at five cm depth was represented by a semicircle with three cm radius. The flux from the source was 10.61 cm h⁻¹ (Eq. (4)), equivalent to a dripper output of 1,200 cm³ h⁻¹. (4) $$\text{σ}=\frac{Q}{4\text{π}{r}^2},$$ where σ is the flux from the surface of the source (cm h⁻¹), Q is the total volumetric flux (cm³ h⁻¹), and r is the radius of the spherical source (cm).

The finite element mesh of the flow domain, which determines the level of model resolution in the calculations, was set to 0.5 cm in both z and h directions. No flux was allowed through the column boundaries. The infiltration source was set as a variable flux boundary so that water and solute applications could be controlled according to the two irrigation scenarios described for the wedge columns (Fig. 2). Nitrate absorption was predicted in the same way as for the horizontal columns. The time-variable boundary condition was used to apply the solute (110 µmol_c NO₃^– cm⁻³) only for the first 0.5 h of water application to the column.

Statistical analysis

The root mean square error (RMSE) was calculated as the error between the measured and simulated water content and NO₃^– concentrations. Comparisons of the RMSE values with predictions from different parameter sets allowed those that produced the lowest errors to be identified. Comparisons of simulated RMSE values with those calculated from measured data allowed the significance of the model error to be assessed in relation to measurement error. This method has been commonly used to measure the quality of model predictions in previous studies (Skaggs et al., 2004; Ajdary et al., 2007; Arbat et al., 2008; Patel & Rajput, 2008). Part of the inverse modelling in Hydrus 1D allows calculation of a correlation matrix that specifies the correlation between the fitted coefficients and statistical information about the fitted parameters.

Parameter determination and model validation in horizontal columns

Values for soil hydraulic parameters estimated from column Set A are presented in Table 2. The two inverse scenarios produced slightly differing values, with the Fit All parameters having a lower value for the objective function Φ than the Set Measured suggesting that the former gave a slightly closer fit between the predicted and measured soil water profiles. The 95% CIs for the optimised values for α, n, and l where smaller compared to the values when all parameters were optimised. The parameters in the Fit All scenario had higher uncertainty and greater correlation between fitted parameters compared to the Set Measured parameters (Tables 2 and 3). The correlation matrix shows there were three high values in the Fit All parameter function compared to one in the Set Measured results; α and n being highly correlated (Table 3, bold entries). High correlation values (magnitude >0.9) indicate parameter nonuniqueness and a correspondingly high uncertainty (Hopmans et al., 2002; Šimůnek & Van Genuchten, 1996). The optimised value of θ_r was 0.112 and the 95% CI that ranged from −0.135 to 0.359, which includes the measured value (0.054 ± 0.003). Saturated water content (θ_s) was 0.56 (0.553–0.567) and was significantly different from the independently measured values (0.58 ± 0.04; Table 2). The measured K_sat (0.104 ± 0.0299) fall within the 95% CI [0.054–0.176; mean best fit value of 0.115] of the optimised value.

The two parameter sets were tested against independently measured data from the horizontal columns (Column data Set B) and the point-source wedge experiments. Predicted water distributions and measured water content in the horizontal columns from column Set B are shown in Fig. 3A. We have plotted θ and NO₃^– profiles against the Boltzmann variable X (distance/√time; cm s^−1/2; Smiles et al., 1978; Phillips & Bond, 1989). Both parameter sets showed similarly good correspondence to the measured data (continuous and dotted lines), which is confirmed by the RMSE and R² values (Table 4). Accurate predictions of θ profiles after absorbing water for 80 and 320 min are not surprising given that data from Set A (27–70 min) are not statistically different from Set B when normalised against the Boltzmann variable. The results, however, do show that under the unsaturated absorption scenario, predictions made for longer times (80 and 320 min) are still accurate when the predictions are extended beyond the range of optimised data. The consistency between predictions based on short time θ and NO₃^– data (Set A) and longer time θ and NO₃^– data (Set B) further demonstrate the experiments had a good degree of repeatability and produced consistent data across a range of time scales.

In comparison, predicted water content using soil hydraulic parameters determined by Rosetta are shown in Fig. 4. These data show the piston front to be significantly behind the measured data (combined data from column Set A and B) especially when the water content at −33 and −1,500 kPa are included in the model. Removal of the FeO before determining the sand, silt, and clay content did not improve the predictions of water retention parameters or saturated hydraulic conductivity. Deriving the water retention parameters and saturated hydraulic conductivity with Rosetta, values that use PTFs to predict the Van Genuchten parameters, produced significant inaccuracies in the predictions (see Fig. 4). In these experiments they are certainly less accurate than parameters derive from inverse modelling (Fig. 3).

The predicted NO₃^– distribution is shown in Fig. 3B. The measured and predicted NO₃^– distributions were compared from both Set A (27–70 min) and Set B (80 and 320 min) because the NO₃^– distribution was not used in the inverse optimisation. The good prediction of NO₃^– distribution by the two parameter sets, when the Langmuir isotherm parameters were included in the simulations, is confirmed by the R² and RMSE values (Table 4). Due to the inaccuracy of water content prediction using Rosetta, no NO₃^– data are presented.

Model validation using 3D wedge infiltration

The predicted distributions of water and NO₃^– throughout the soil wedge after the two irrigation scenarios are shown in Fig. 5. This figure also shows the positions of horizontal and vertical transects presented in Figs. 6 and 7. These figures show the agreement between the measured and predicted water and NO₃^– profiles in the wedge. Comparisons of the RMSE and R² calculations indicated that both the Fit All and Set Measured parameter sets predicted very similar distributions, although the Fit All parameters produced slightly better predictions of NO₃^– distribution in the longer irrigation scenario. The measured RMSE, calculated from columns where duplicate measurements were taken at identical times and locations, were similar to the RMSE of the predicted values (Table 5). This indicates that the errors between the measured and predicted values were very similar to the errors of replicate measurements at identical points and times in the wedge experiments. The predictions achieved using the two parameter sets were therefore considered to be suitable to estimate water and NO₃^– distribution in the point-source flow scenario of the wedge column. The “Set Measured” parameters are preferred because there is less auto correlation between the fitted parameters (α, n, l).