Optimal Choice of Soil Hydraulic Parameters for Simulating the Unsaturated Flow: a Case Study on the Island of Miyakojima, Japan

We examined the influence of input soil hydraulic parameters on HYDRUS-1D simulations of evapotranspiration and volumetric water contents (VWCs) in the unsaturated zone of a sugarcane field on the island of Miyakojima, Japan. We first optimized the parameters for root water uptake and examined the influence of soil hydraulic parameters (water retention curve and hydraulic conductivity) on simulations of evapotranspiration. We then compared VWCs simulated using measured soil hydraulic parameters with those using pedotransfer estimates obtained with the ROSETTA software package. Our results confirm that it is important to always use soil hydraulic parameters based on measured data, if available, when simulating evapotranspiration and unsaturated water flow processes, rather than pedotransfer functions.


Introduction
There are no rivers on Miyakojima, a semitropical island in southernmost Japan, because of its flat topography and the high permeability of the limestone that forms the island.Local residents there are dependent on groundwater for almost all of their domestic water use.In the Okinawa region, temperatures are predicted to rise in response to climate change, while annual rainfall is expected to decrease [1], with resultant depletion of groundwater resources becoming a concern.The main land use on Miyakojima is sugarcane farming; thus, understanding both water movement in the unsaturated zone of the farmland soil and the total water budget is important.
The HYDRUS-1D software package [2] has often been used for analyses of these type of problems [3,4].HYDRUS-1D provides versatile numerical modeling of the movement of moisture, solutes, and heat in soil.One option in the code is to estimate soil hydraulic properties by using pedo-transfer functions (PTFs).Since it is difficult to measure soil hydraulic parameters, PTFs that estimate them from readily measurable soil characteristics, such as particle size distribution and bulk density provide a very attractive tool for numerical analyses.
The accuracy of HYDRUS-1D simulations has been analytically verified [5].Most studies to evaluate the performance of HYDRUS-1D were simulations of the transfer of heat and moisture in semiarid and humid regions.For example, Saito et al. [6] reported that HYDRUS-1D was useful for predicting the transfer of heat and moisture in sandy soils, and Kato et al. [7] reported that it was useful for predicting soil temperature and moisture in volcanic soils.
Other studies have used vadose zone models and PTFs.Steinzer et al. [8] simulated evapotranspiration, infiltration, and VWCs distribution in lysimeter experiments on sandy and clay soils in Germany.Wang et al. [9] simulated groundwater recharge in sand and sandy loam soil by using HYDRUS-1D.They found that recharge in their soils was strongly dependent on the parameter n of the van Genuchten model [10] and uncertainties in the simulated recharge were affected by uncertainties in the n estimated from PTFs.
Thus, although HYDRUS-1D is known to be useful for simulations of water movement in the unsaturated zone, the code must be tested in a particular area before practical applications; for example, for the development of a water management plan.Wang et al. [11] simulated groundwater recharge at four sites in the continental United States with different climate conditions using HYDRUS-1D along with datasets for sand and loamy sand.They showed that the distribution patterns of mean annual groundwater recharge varied considerably across the sites, mainly depending on soil texture and climatic conditions.
To date, the use of HYDRUS-1D in the island of Miyakojima has not yet been tested.Although the necessary weather data are readily available from the Japan Meteorological Agency, the collection of soil data is more difficult.Most of the soil data collected and analyzed in past investigations for land development projects consists only of particle size distributions and bulk density; water retention curves were not always determined.Accurate simulation results are dependent on the quality of the soil hydraulic parameters used as input to the simulation.Okamoto et al. [12] reported that the retention curve of Shimajiri mahji soil (dark-red soil [13], which was classified as a Cambisol [14]) estimated using the ROSETTA software package [15][16][17] was considerably different from that derived from measured data.However, they did not comment on the influence of this difference on the simulated water budget.
To validate the use of HYDRUS-1D on Miyakojima, we examined the influence of input soil hydraulic parameters on HYDRUS-1D simulations of evapotranspiration and VWCs.We first optimized the parameters for root water uptake and examined the influence of soil parameters on simulations of evapotranspiration (details are shown in the later section).We then compared VWCs simulated using measured soil hydraulic parameters with those simulated using parameters derived by using ROSETTA software.
We measured soil particle size distributions (%), bulk densities (g•cm −3 ), hydraulic conductivities (cm•day −1 ), and water retention curves in the laboratory.Soil particle size distributions were obtained using sieve analysis for particle sizes greater than 75 μm and by hydrometer analysis for particle sizes smaller than 75 μm, according to Japanese Industrial Standards (JIS A1202).We measured bulk densities by using the 100 cm 3 soil core sampler, hydraulic conductivities by the constant-head method, and water retention curves by the three methods shown in Table 1.
Table 1.Experimental procedures used for VWC measurements at selected soil suction ranges.

Overview of HYDRUS-1D
One-dimensional water flow in soil is described by the Richards equation as follows: where θ is volumetric water content (cm 3 •cm −3 ), t is time (s), z is the spatial coordinate, assumed positive upward (cm), h is pressure head (cm), S is a sink term for water uptake by plant roots (cm 3 •cm −3 •s −1 ), and K is unsaturated hydraulic conductivity (cm•s −1 ).We applied the van Genuchten-Mualem equation [10,20] (VG hereafter) to estimate the water retention curve and hydraulic conductivity of the soil by using Equations ( 2) and (3), respectively: ( ) ( ) where θ r is residual water content (cm 3 •cm −3 ), θ s is saturated water content (cm 3 •cm −3 ), and S e is normalized water content.α (cm −1 ), n, m (= 1 − 1/n) and l (= 0.5) are empirical parameters, and K s is saturated hydraulic conductivity (cm s −1 ).
We simulated five soil layers (0-20, 20-40, 40-60, 60-80, and 80-120 cm) that correspond to the depth ranges of the soil samples we collected.Atmospheric conditions of daily precipitation and potential evaporation were used as upper boundary conditions.Potential transpiration was used to calculate water uptake by plant roots.The lower boundary condition was free drainage.Initial VWC was the VWC measured in the field.

Crop Model in HYDRUS-1D
We used the method of van Genuchten et al. [21] to simulate water uptake by plant roots: ( ) where α(h) is the water stress response function, S p is potential water uptake by plants roots, h 50 is the soil suction at which water uptake by roots is reduced by 50%, and p is an empirical component that is usually assumed to be 3 [22].We used p = 3 to optimize h 50 .
To determine daily root length, we used a logistic root growth function in HYDRUS-1D [2] and the results of a previous study [23].We used the leaf area index (LAI) growth model of Larsbo and Jarvis [24].The LAI growth model from crop emergence day (D min ) to the day of maximum LAI (D max ) is described by Equation ( 6), and that from D max to the day of harvest (D harv ) by Equation ( 7): ( ) ( ) where LAI min is LAI at D min , LAI max is LAI at D max , LAI harv is LAI at D harv , D* is the number of days after planting, and x 1 and x 2 are empirical components.We derived daily values of LAI from the results of a previous study of sugarcane [25] (Figure 2).Potential evapotranspiration (ET p ) was calculated using the Penman equation ( 8) [26,27]: ( ) where ET p is the reference evapotranspiration (mm•day −1 ), Δ is the slope of the saturation vapor pressure curve (kPa•°C −1 ), γ is the psychometric constant (kPa•°C −1 ), S is the net radiation (MJ•m −2 •day −1 ), L is the latent heat of evaporation (MJ•kg −1 ) is, u 2 is wind speed at 2 m height (m•s −1 ), (e sa − e a ) is the saturation vapor pressure deficit (kPa).ET p was partitioned into potential evaporation (E p ) and potential transpiration (T p ) using Campbell's equation [28]: Total precipitation during the observation period was 675.5 mm (Figure 3).Since the days after planting during the observation period varied from 189 to 313 and LAI from 2.4 to 3.3 (from Equations ( 6) and ( 7)), E p was low throughout the observation period.

Comparison of Measured Soil Hydraulic Parameters with Those Estimated Using ROSETTA
We determined the parameters for the VG equation by using our analyses of soil samples and by application of the ROSETTA module in HYDRUS-1D.We determined VG parameters so as to minimize the difference between measured values and estimated values according to the nonlinear least-squares method using the solver of an EXCEL add-in.ROSETTA can use combinations of soil texture, particle size distribution, bulk density, and one or two points on the water retention curve to determine the parameters for the VG equation and for hydraulic conductivity.In this study, we used particle size distribution and bulk density, since these properties were readily available.
Measured particle size distribution, bulk density and hydraulic conductivity are shown in Table 2.For all samples, the sand content was less than 20% and the soil texture was silt.There were no clear differences in bulk density among soil layers.The standard deviation of hydraulic conductivity was large for all layers, which we attributed to cracks formed in response to shrinkage during drying of the soil samples.We, therefore, used the geometric mean of hydraulic conductivity in the simulations.In our application of HYDRUS-1D, we considered four combinations of input soil hydraulic parameters (Cases 1 to 4; Table 3).

Simulation of Evapotranspiration and Volumetric Water Contents
To optimize the value of h 50 , we needed to minimize the difference between measured and simulated evapotranspiration for each of Cases 1 to 4. To achieve this, for each case we ran HYDRUS-1D using h 50 values from 100 to 1000 cm at intervals of 100 cm.Then, we chose the optimum value of h 50 as the value with the lowest root-mean-square error (RMSE) for total evapotranspiration calculated at 10-day intervals.
To examine the influence of input soil hydraulic parameters on simulated VWCs movement, we calculated the RMSE between measured VWCs and simulated VWCs (optimized h 50 ) for each of the four soil layers for Cases 1 to 4.

Comparison of Measured Soil Hydraulic Parameters with Those Estimated by ROSETTA
In general, ROSETTA underestimated measured VWCs, particularly for soil suctions greater than 1000 cm (Figure 4).Okamoto et al. [12] reported similar results from their application of ROSETTA to Shimajiri mahji soil.Since Jahgaru soil has a poorly developed structure [19] it often does not drain well [18].Therefore, the measured VWC values at high suction tended to be larger than that estimated by ROSETTA (Figure 4).Schaap and Leij [29] reported that the use of PTFs to estimate soil hydraulic parameters might depend strongly on the data used for calibration.The characteristics of Jahgaru soil might be different from the soil data used to develop on the parameters used in ROSETTA.4).We attributed the difference in these results to the influence of cracks and aggregations of soil in the sugarcane field.Note: ROS = from ROSETTA.

Optimization of h 50 and Simulation of Evapotranspiration
The simulations run to optimize h 50 (Figure 5) show that total evapotranspiration (TET) increased with increasing h 50 for all four cases considered.The optimal h 50 (smallest RMSE) was 600 cm for Case 1, 400 cm for Case 2, 300 cm for Case 3 and 700 cm for Case 4. For each of the four cases, the TET estimated with the optimized h 50 was almost the same as the measured TET (204.5 mm).The simulation results for the pairs of cases with the same retention curve (Cases 1 and 4, Cases 2 and 3) were similar (compare Figures 5 and 6).The suction required to deplete VWCs during normal growth has been reported to be about 1000 cm [30]; therefore, we considered that Case 1 (600 cm) and Case 4 (700 cm) provided the more realistic values of h 50 for application in HYDRUS-1D, even though the RMSEs of Cases 2 and 3 were lower.The common factor for Cases 1 and 4 was the use of the measured retention curve.
Thus, comparison of our simulation results for evapotranspiration did not clearly indicate which soil hydraulic parameters were better for application in HYDRUS-1D.

Influence of Soil Hydraulic Parameters on Volumetric Water Contents Simulation
The time series of simulated VWCs (Figure 7) and RMSEs between measured and simulated VWCs (Table 5) show that the simulated VWC values for Cases 2 and 3 (retention parameters estimated using ROSETTA) were lower than measured values, whereas for Cases 1 and 4 (retention parameters calculated from measured data) simulated and measured VWCs agreed well.When considering the water budget, water recharge (Re) is calculated using where P is precipitation, ET is evapotranspiration and ΔSW is the change of VWCs.To estimate water recharge from this equation, accurate simulations of VWCs and evapotranspiration are needed.The results of our simulations of evapotranspiration for Cases 1 to 4 (Figure 6) did not differ greatly (Section 3.2.1).However, our simulations of VWCs indicated that Cases 1 and 4 provided better results than Cases 2 and 3.These results indicate the importance of using soil hydraulic parameters based on retention curves derived from measured data to calculate the water budget for the island of Miyakojima, and likely for many or most other applications.

Conclusions
We drew the following main conclusions about the optimum application of HYDRUS-1D in our study area.
2) Optimized values of h 50 were dependent on the parameters defined by the retention curve.
Simulated and measured total evapotranspiration rates agreed well for all four cases considered.Since, for normal growth the amount of suction required to deplete VWCs is about 1000 cm, we consider the h 50 values we obtained that were closest to 1000 cm to be the more realistic.Thus, our HYDRUS-1D simulations using the measured soil hydraulic parameters provided better results than those based on parameters estimated by ROSETTA.3) VWCs simulated by HYDRUS-1D using parameters estimated by ROSETTA were lower than the measured values, whereas those using measured parameters agreed well with measured values.
Our study confirmed that it is important to use soil hydraulic parameters derived from measured retention data on the island of Miyakojima, rather than estimates obtained with pedotransfer function.

Figure 6 .
Figure 6.Comparison of simulated, measured and potential evapotranspiration at ten-day intervals (MET) from 29 August to 31 December 2009.

min. Measurement started on 29 August 2009 and ended on 31 December 2009. Cultivation started on 21 February 2009 and ended on 20 January 2010.
Groundwater basins on the island of Miyakojima.

Table 2 .
Measured particle distribution, bulk density, and hydraulic conductivity of all layers.

Table 4 .
Parameters used in simulations.