Applicability of Hydrus-1D in a Mediterranean Mountain Area Submitted to Land Use Changes

The aim of this research is to evaluate the reliability and accuracy of the Hydrus-1D model to simulate the measured dynamics of water ﬂow in a silt loam soil proﬁle located in an abandoned crop area. The paper includes a physical and chemical characterization of the soil, and hydraulic properties characteristics as well. Several techniques and devices were used to develop the experiment in both, ﬁeld and laboratory scales. The last part of the study was the Hydrus-1D simulation using real rain events and evapotranspi-ration rates. In summary, it could predict accurately the water dynamics of this “natural” scenario.


Introduction
During the last decades, numerical models have been increasingly used to predict and to analyse water flow and solute transport in the unsaturated zone. Although a number of analytical approaches have been developed to solve water flow equations [1,2], conceptual and numerical difficulties still exist, especially at the transient scale flow and with multidimensional field applications [3].
Soil macropores have been studied as a very important factor in the movement of water in both vertical and lateral directions [4][5][6]. Soil porosity largely determines soil hydraulic conductivity [7], usually represented by a lognormal function related with a high heterogeneity of pore size, and showing high increases when water content into the porous media increases.
Quantitative studies of soil water transport due to macropore flow have resulted in several empiric and theoretical approaches aimed at describing such process. Physical models, usually derived from Darcy's law, for example, [8,9], including some of these as are Green and Ampt [10], Philip [11], and van Genuchten [12], which explain different patterns of water flow into the soil matrix.
Following authors such as Skopp [13], Flühler et al. [14] or Jarvis [15] can observe that water flow into the soil presents certain instability, which can be due to the characteristics of the flow itself, and to structural features of the soil. The latter include water repellence, the slope of the different layers in the soil matrix [16], and swelling-shrinking processes, which involve a nonrigid soil as well.
Some of these unsteady water flows, mainly the preferential flow, can be especially important in soils with high swelling-clay content, enhancing the risk of pollution by agricultural chemicals when the hydraulic conductivity values increase. In addition, the specific characteristics of the local area also influence water pathways. Christiansen et al. [17], in a recent study, showed that preferential flow varied significantly with topography and groundwater within a basin.
Many field studies have shown the importance of preferential flow in crop fields and areas with artificial drainage networks. This flow is described by different conceptual models based on the classic Richards equation [9], which explain the behaviour of a vertical transient flux in an unsaturated area. The most useful models are double porosity and multi-porosity and double permeability and multipermeability [15,[18][19][20]. These models assume that the porous media is divided in two different regions, one associated at the interaggregate porosity and the other at the intraaggregate porosity. Hence, we would like to emphasize some of these models, for example, SWATRER [21], SWMII [22], CHAIN-2D [23], and especially HYDRUS-1D [24], which will be used in this study.
Land use changes involve landscape and vegetation shifts, which ultimately influence the physical characteristics of hillslopes and soils [25,26]. Tillage techniques involve soil alterations which can lead to an increased water retention capacity [27,28]. The presence of certain plant species can enhance soil organic matter hydrophobic properties [29]. Structural stability may also vary with land use changes, involving a degradation of this soil property, and decreasing the permeability and porous media volume [30,31]. The abandonment of cultivation in montane areas usually results in the increase of meadow and forest surface. The changes in the vegetation suggest changes in the organic matter content, structural stability, and aggregate size, and hence in bulk density and porosity, especially at surface soil level according to [32,33]. Also, these changes have shown an increase in soil microbial activity, which influences soil hydrodynamic properties such as infiltration rates and hydraulic conductivity [34,35].
Afforestation processes following agricultural abandonment in mountain areas due to socioeconomic changes have led to a significant increase in forest cover throughout upland areas in Northern Spain. Geomorphological and hydrological processes associated to these land-use changes have been intensively studied at the Vallcebre experimental basins (Eastern Pyrenees, NE Spain). At this site, the terraced topography influences the spatial patterns of the soil moisture [36] and runoff processes [37][38][39]. Besides, a 25% increase in forest cover during the second half of the XXth century [36,40] may indirectly affect soil physical properties and therefore water movement in the soil.
The objective of this study is to evaluate the reliability and accuracy of the Hydrus-1D model [24] to simulate the measured dynamics of water flow in a silt loam soil profile located in a terraced hillslope whose cultivation was abandoned during the XXth century and which was overgrown with pasture. For this purpose, the main objective has been addressed in four distinct tasks: (i) to determine the hydraulic properties of the soil profile, (ii) to parameterise the van Genuchten model and field saturated hydraulic conductivity, (iii) to calibrate Hydrus-1D using the matric potentials obtained for the initial boundary conditions; (iv) to evaluate Hydrus-1D model using water contents and pressure heads.

Materials and Methods
The Can Vila research basin (0.56 km 2 ) is located in the head basin of the Llobregat River, northeast of Spain (42 • 12 N; 1 • 49 E) ( Figure 1). The topography of Can Vila basin presents a modified hill slopes series in terraced landscape, with an altitude between 1100 and 1600 m a.s.l. The climate is warm temperate (Papadakis classification) with a mean annual temperature of 7.3 • C and an average annual precipitation rate of 924 mm, being spring and autumn the wettest seasons. The evaporation and transpiration rate is 700 mm  [41]. The climactic vegetation of the study area is Quercus pubetcens forest although at present Aphyllantion-type calciophilous grassland [42] is dominant in the terraced areas, together with spontaneous afforestation patches of Pinus sylvestris with a poor Buxus sempervirens shrub cover [43]. The soils, which can be classified as Ochrepts, Orthents, and Aquepts according to soil taxonomy system [44], present a high spatial variability, especially in the infiltration rates. During wetting-up conditions, saturation patches appear in the inner part of the terraces, underlain by mudstones with low-infiltration capacity which appeared after terrace construction. These terraces are drained by an artificial channel network [38,45].

Hydraulic Properties of the Soil Profile.
Our experiment was carried out at La Call terrace. The experimental plot is an old abandoned crop terrace located in Can Vila basin. The experimental dataset included water potential data measured at 0.2, 0.4, and 0.6 meters of depth using SKT600 tensiometers (Skye Instruments Ltd.). The data was recorded every 20 minutes by a data-logger DT-500 (DataTaker). For the observed water contents, two time domain reflectometry (TDR) profiles (A and B) were used. The profile A was located in the inner part of the terrace and profile B in the middle part of the same terrace. Water content data were collected between surface and 0.6 meters of depth, using 3 TDR probes of 0.2 meters long per profile and a Tektronics 1502C cable tester, for measuring the data. Also, a set of unaltered samples (100 cm 3 ) from a soil profile of the experimental plot were used to characterize the physical, chemical, and hydraulic soil properties. The samples were obtained between surface and 0.85 meters of depth.
Volumetric water content of unaltered samples was determined using the "sand box" method for potentials between saturation and −20 kPa [46] and the pressure membrane for potentials between −100 and −1500 kPa [47]. Water retention curves were fitted to the van Genuchten model, using the RETC code [48] obtaining θr, θs (residual and saturation volumetric water content, resp.), α, and N (shape parameters). In addition, the following properties were determined: particle size distribution, bulk density, porosity, calcium carbonate, and organic matter content.
In the same plot, we determined the field-saturated hydraulic conductivity (Kfs), as well. The observed Kfs data were obtained for 3 depths (0.15, 0.25, and 0.50 meters) and 3 replica per depth, using for it the Guelph permeameter method [49]. The diameter of the well was 6 cm. The rainfall event dataset was obtained with an automatic pluviometer with 0.2 mm of accuracy, and the reference evaporation-transpiration was calculated according to Penman-Monteith equation [50]. Root distribution was based on field observations and the trapezoidal uptake function [51]. The soil water stress reaction was calculated using the S-shape function [52]. The expression is as follow where S(h) is the volume of water removed from a unit volume of soil per unit time due to plant water uptake, h φ is the osmotic head, β(h) is a prescribed dimensionless function of the soil water pressure head (0 ≤ β ≤ 1), S p is the potential water uptake rate, h 50 represents the pressure head at which the water extraction rate is reduced by 50% during conditions of negligible osmotic stress, and p is a experimental constant.  50 parameter from the S-shape function was −359 cm, and the p-parameter was equal to 3. In order to improve the simulation of volumetric water contents, the closed form of the van Genuchten model with m = 1 − 2 · n −1 was used. The model is expressed as

Hydrus-1D Calibration and
where ψ is the soil water pressure head, α and N are curve shape parameters, and Se is the relative saturation that is expressed in actual, residual, and saturated volumetric water content (θ, θr, θs) as Field-saturated hydraulic conductivities were the mean geometric values obtained for both seasons (wet and dry) at different sampling levels (Table 3). For the first layer (from surface to 0.8 m depth), the Kfs value was 1.64 m·day −1 , and for the second layer (from 0.8 cm to bottom profile), the Kfs value was 0.16 m·day −1 .
To validate the water content data using Hydrus-1D, the volumetric water contents calculated according to Topp et al. [53] were used. These data were obtained for the two TDR soil profiles. Table 1 shows the soil properties of the studied profile. The textural class of the profile was classified as silt loam (according to USDA), with silt content always higher than 57 g·kg −1 , sand content between 110 g·kg −1 and 210 g·kg −1 , and clay content between 200 g·kg −1 and 280 g·kg −1 . Mean bulk density was 1.38 Mg·m −3 increasing with depth and yielding a value for total porosity about 47%. Mean organic matter content, about 40 g·kg −1 , indicates a high organic carbon content decreasing with depth. Mean calcium carbonate content is reasonably high, about 320 g·kg −1 . According to coefficient of variation, the behaviour on the whole soil profile can be considered as a low variability (see Table 1).  Table 2 presents different values obtained with the van Genuchten equation, note that N-parameter is less than 1.20 except for the deepest sampling level, for which the corresponding value is 1.26. Despite to obtain an excellent fitted values for Nparameter, these values were solved for Hydrus-1D using an  air entry value −2 cm with m-parameter, this fact improved the solution of the problem. Figure 2 shows, at daily scale, a simulation for calibration starting in dry conditions, with an accumulated precipitation of about 264.4 mm after a dry period. There were significant differences between the observations and the values predicted by Hydrus-1D. The simulations revealed a significant time delay at deeper levels (0.4 and 0.6 m depth), immediately after the calculations started. The response of the model was slower than that of the measured data, especially at 0.6 cm depth. The model predicted correctly (or only slightly overestimated) observed pressure heads near saturation, but was unable to predict pressure heads during periods without rainfall when water uptake was important, due to high evaporative demand. Pressure head at 0.2 m depth fitted well throughout the simulated period, showing only a slight underestimation at the start of the simulation, after the first rainfall event. The two tensiometric observation points (0.4 and 0.6 m depth) showed a time delay at the start of the simulation. At this moment, the model did not describe well the soil's fast response immediately after wettingup. On the other hand, Figure 3 shows a simulation for the same period starting with wet conditions. The results indicated that the response of the model is faster than that measured after initial dry conditions, especially at the deeper levels. An acceptable simulation of the soil wetting-up and transition to saturation was observed, improving the results obtained using dry initial conditions and decreasing the time delay between the first rainfall event and the model response, even though the drying out of the curve is not well represented by the model.

Simulating Water Contents with Hydrus-1D.
The results have shown that Hydrus-1D predicted reasonably well volumetric water content using an air entry value of −2 cm in the van Genuchten equation, obtaining correlation coefficients of about 0.75 (P ≤ 0.05). The relationship between measured and predicted mean volumetric water contents data for both TDR profiles used for calibration (Figure 4) showed that Hydrus-1D gave an acceptable fit from surface to 0.6 meters of depth. The TDR profile in the inner part of the terrace (profile A) showed the most accurate fit, with less dispersion and a similar trend throughout the simulated period. Profile B located in the outer part of the terrace, and the water contents were overestimated, with larger overestimations after rainfall events about 50 mm·day −1 . However, the model described acceptably the wettingup transition period after the first rainfall event.

Discussion and Conclusions
The studied soil profiles showed high water content values (at saturation ∼ =0. 50  This feature may be related to the high silt and clay content, greater than 75 g·kg −1 for all soil profile [54][55][56][57]. Organic matter content close to surface was about 7 g·kg −1 [55,58,59]. The information obtained by several soil samples (3 replica per level) of the seven sampling levels for a silt loam soil profile did not improve the Hydrus-1D simulation, because it only used two different materials (or layers) related with the soil textural class. The values of saturated hydraulic conductivity found in the field test could be a consequence of macropores and preferential flows; note that these soils have high swelling-clay content. In addition, high calcium carbonate content gives a sustaining structure during the macropores formation [60]. This fact would explain the increase in the water contents in the deepest layers, being preferential flow controlled by a macroporous network [61]. Swelling possibly closed the macropores causing a fast reaction of the soil which the model Hydrus-1D could not describe well, needing more time to acquire stability during the simulation. This could explain the poor adjustment of the model during the initial time steps of prediction after dry initial conditions. However, when the soil profile was already moist, the model response improved, not showing the time delay of the deeper levels.
In summary, the simulation of Hydrus-1D using the van Genuchten equation with an air entry value of −2 cm was adequate to estimate soil water contents and pressure heads for a silt loam soil profile located in Can Vila basin, in spite of the differences obtained between observed and predicted data. Our results are comparable to those obtained by other authors [62,63]. To conclude, we wanted to ascertain that the algorithm of Hydrus-1D solved correctly the Richards equation for this silt loam soil profile under natural conditions. Although the model simulated the pressure head data with some differences, especially in the dryout events, it should be noted that Hydrus-1D is a robust model to predict soil water contents on different scenarios with transient conditions.