Abstract
A FORTRAN code for liquid water flow in unsaturated soil under the isothermal condition was developed to simulate water infiltration into Yolo light clay. The governing equation, that is, Richards’ equation, was approximated by the finite-difference method. A normalized sensitivity coefficient was used in the sensitivity analysis of Richards’ equation. Normalized sensitivity coefficient was calculated using one-at-a-time (OAT) method and elementary effects (EE) method based on hydraulic functions for matric suction and hydraulic conductivity. Results from EE method provided additional insight into model input parameters, such as input parameter linearity and oscillating sign effect. Boundary volumetric water content ( (upper bound)) and saturated volumetric water content () were consistently found to be the most sensitive parameters corresponding to positive and negative relations, as given by the hydraulic functions. In addition, although initial volumetric water content ( (initial cond)) and time-step size (Δt), respectively, possessed a great amount of sensitivity coefficient and uncertainty value, they did not exhibit significant influence on model output as demonstrated by spatial discretization size (Δz). The input multiplication of parameters sensitivity coefficient and uncertainty value was found to affect the outcome of model simulation, in which parameter with the highest value was found to be Δz.
1. Introduction
Sensitivity analysis is used for various reasons, such as decision-making or development of recommendations, communication, increasing understanding or quantification of system, and model development. In model development, it can be used for the purposes of model validation or accuracy, simplification, calibration, and coping with poor or missing data and even to identify important parameter for further studies [1].
More than a dozen sensitivity analysis methods are available, ranging from one-at-a-time (OAT) to variance-based methods [2, 3]. In a fundamental level, sensitivity analysis is a tool to assess the effect of changes in input parameter value on output value of a simulation model. In this aspect, the sensitivity coefficient, in a normalized form, is given in the following relation: where is referred to as normalized sensitivity coefficient for th input parameter at th observation point, is model dependent variable value at th observation point, and is th input parameter value. This method utilizes derivative at a single point and similarly it can be applied as OAT method when one input parameter is varied while holding other parameters fixed. However, the former and the latter methods do not explore other input space factors in which more than one input parameter is varied. Despite this disadvantage, Saltelli and Annoni [4] noticed that researchers continuously practice OAT method mainly due to few advantages claimed, for example, a safe starting point where the model properties are well known and all OAT sensitivities relative to a starting point. Although variance-based method is the best practice, Saltelli and Annoni [4] have suggested the use of elementary effects method, which is an enhancement of OAT method, when computation time is expensive, for instance, in numerical simulation that is computationally demanding.
Elementary effects method is accomplished through the use of a technical scheme to generate trajectories. Each trajectory consists of a number of steps in which each step is referred to an increment or decrement of an input parameter value. The base condition for each trajectory is different from the others, and it is selected randomly. The random version of trajectory generation is as follows [5]:where is generated trajectory in the form of matrix with dimension , where is the number of independent input parameters, is a value in and is the number of levels, is matrix of 1’s, is a randomly chosen base value, is lower triangular matrix of 1’s, is -dimensional diagonal matrix in which each element is either +1 or −1, by random generation, and is -by- random permutation matrix that each row and column of the matrix with only one element equal to 1 while the other elements of the matrix are zero.
The generated trajectories can be screened to obtain a subset of trajectories with the greatest geometric distances. The trajectories scanning to maximize geometric distances between all the pairs of points between two trajectories is as follows [6]:where is distance between a pair of trajectories, and , is th coordinate of the th point of the th trajectory, and is th coordinate of the th point of the th trajectory.
The sensitivity coefficient of an input parameter in elementary effects method is presented as , which is the mean of elementary effects (). is the mean of absolute values of the elementary effects, which is used to avoid cancellation of difference signs in the mean value. The sensitivity measures (, , and ) and are given by [5]:where and are simulation result before and after increment or decrement of value, that is, , which can be either of positive or negative value, is the total number of trajectories, is elementary effects of input parameter at trajectory, and is standard deviation of input parameter.
The aim of the current work is to carry out sensitivity analysis on water infiltration into unsaturated soil as governed by Richards’ equation, that is, governing equation of soil water flow, and use it as an evaluating method to validate the simulation source code with analytical solution. Thus, the objectives of this study are to (1) determine the sensitivity coefficient and (2) to validate model simulation with Philip’s semianalytical solution from literatures using the sensitivity coefficient, under a hypothetical assumption. In this study, we used the water infiltration results from Haverkamp et al. [7] and Kabala and Milly [8] to verify the simulation.
2. Materials and Methods
2.1. The Governing Equation of Water Flow in Unsaturated Soil, and Its Numerical Solution
The governing equation for transient liquid water flow in soil may be described as [9]:where is volumetric water content (m3 m−3), is time (s), indicates vertical distance (m), is hydraulic conductivity of soil (m s−1), is matric pressure head (m), is vector unit with a value of positive one when it is vertically downwards.
Equation (8) was approximated numerically and its algebra was implemented in FORTRAN 2008 using Simply FORTRAN Integrated Development Environment. The spatial discretization method used is termed as cell-centered finite-difference, and the temporal discretization method used was the fully implicit scheme. In order to avoid unnecessary redundancy, we only provide the algebra for (8) that is used for sensitivity analysis in the current study as follows:where indicates a cell-centered number in -direction in Cartesian coordinate system, (s) is time-step size, (m3 m−3) and (m3 m−3) are volumetric water content at old time level and new time level , respectively, (m s−1) is hydraulic conductivity at the interface between cells and , (m s−1) is hydraulic conductivity at the interface between cells and , is partial derivative of with respect to at the interface between cells and , is partial derivative of with respect to at the interface between cells and , (m), (m), and (m) are corresponding to the spatial sizes of spacing of cells , , and , respectively, (m3 m−3), (m3 m−3), and (m3 m−3) are the volumetric water contents at new time level of cells , , and , respectively. Equation (8) was numerically solved by a fully implicit cell-centered finite-difference scheme without any linearization. An iterative method was used to solve the mathematical algebra of (9), that is, Jacobi iteration [10]. For comparison purpose, modified Newton-Raphson method was also implemented [11]. A convergence factor criterion was used to indicate the condition for iteration termination process, that is, absolute maximum difference for every single cell.
2.2. The Constitutive Functions of Matric Pressure Head () and Hydraulic Conductivity ()
The hydraulic functions used were adopted from Haverkamp et al. [7]:where , , , and are fitting parameters, (m3 m−3) is residual volumetric water content, (m3 m−3) is saturated volumetric water content, and (m s−1) is saturated hydraulic conductivity.
2.3. Numerical Experiment and the Default Setting of Input Parameters of the Flow Problem
Water infiltration into Yolo light clay was used in the numerical experiment. The hydraulic functions for the soil (see (10)), and the coefficients values are shown in Table 1. Initial condition for the volumetric water content was 0.2376 m3 m−3. Lower boundary was set as free-drainage to water flow. Upper boundary was set at 0.495 m3 m−3. After considering the mass balance ratio [9] and iteration number, the time-step size, spatial discretization size, and convergent value were set at corresponding values of 500 s, 1 cm, and 10−12 m3 m−3, respectively. The iteration methods of Jacobi and modified Newton-Raphson were compared. It was found that the minimum iteration number from the latter was equivalent to the iteration number from the former, when the relaxation factor of the latter was set to unity (data not shown). Reducing the relaxation factor from unity would result in increasing iteration number. The numerical solution of (9) did not exhibit convergent problem; thus, Jacobi iteration method was sufficient.
2.4. Statistical Measures
In order to determine the goodness of fit between reference data and simulated results, one statistical equation was implemented. The equation is called absolute residual errors (MA) as follows [12]:where is the simulated data at cell and is the analytical solution as reference data at cell .
3. Results and Discussion
3.1. Simulation Results and Their Accuracy
Based on the conditions as stated in previous section, water infiltration into Yolo light clay was simulated up to 105 s. Data on Philip’s semianalytical solution were collected from Haverkamp et al. [7], hereafter referred to as Philip(H). Simulation results were compared with the data to verify the simulation (Figure 1). It was evident that the simulation results slightly underpredicted the infiltration front of water flow.
To further reinforce the previous claim, some data were extracted from Kabala and Milly [8], as indicated by Philip(K) as in Figure 1, for further comparison. Figure 1 shows that there was a small difference between Philip(K) and Philip(H), but the former was relatively closer to the simulation results than the latter. At this point of observation, we were unable to determine which of the solutions, that is, Philip(K) and Philip(H), provided from the literature was accurate. However, results from the figure clearly indicate that the simulated result was lesser than Philip’s semianalytical solution. Therefore, sensitivity analysis was carried out to determine the sensitivity coefficient for all input parameters and use the sensitivity analysis results to assess the model simulation based on the assumption that possibly the cumulative effect of input parameters, in terms of significant digits approximation, could be contributing to the underprediction of the volumetric water content of the simulation. In addition, sensitivity analysis is one of the most important steps in evaluating the effect of input parameter on simulation results, and it is also used by others for model validation [13–16].
3.2. Sensitivity Analysis and Simulation Model Validation
Negligible sensitivity response could be due to too small perturbation size, and inaccuracy in sensitivity response could be due to too large perturbation size [17]. Values of input parameters were subjected to a perturbation size between −5% and 5% as suggested by Zheng and Bennett [12], and in considering the simulation time, we limit the sensitivity analysis to a simulation time of 105 s. The sensitivity analysis study was based on a single perturbation size of increment or decrement in each simulation. The sensitivity analysis was carried out based on the hydraulic functions, (10), from Haverkamp et al. [7].
There were seven input parameters from Haverkamp hydraulic functions as listed in Table 1. Additional four input parameters were also tested, that is, initial volumetric water content ( (initial cond)), boundary volumetric water content (upper cond), time-step size , and spatial spacing size . The depth at 15.5 cm from the ground surface was used for observation.
The normalized sensitivity coefficients are shown in Figure 2. Generally, there are two groups of sensitivity coefficients, that is, positive and negative relations. In positive relation group, the boundary volumetric water content had the highest sensitivity coefficient. This was followed by initial volumetric water content and saturated hydraulic conductivity. The smallest sensitivity coefficient in the group was the residual volumetric water content. In negative relation group, saturated volumetric water content had the highest sensitivity coefficient, and this group ended with spatial spacing size and time-step size as the smallest sensitivity coefficient.
For comparison purpose, elementary effects method was also used to calculate normalized sensitivity coefficient. We assumed only random generation in -dimensional diagonal matrix (), and then (2) was used to generate 50 trajectories. Equation (3) was used to screen out 4 trajectories with the greatest geometric distance of those trajectories. Equations (4) to (7) were used to calculate the elementary effects, mean of elementary effects, mean of absolute values of the elementary effects, and standard deviation, respectively. The mean of elementary effects was modified to calculate the normalized sensitivity coefficient. The results are shown in Table 2. The sensitivity coefficient has identical ranking as those obtained in Figure 2, except for the coefficient of input parameter. Similar values of and indicate linear effect on few input parameters in positive (, , and (initial cond)) and negative ( and ) relations. Other input parameters have shown the effect of oscillating sign that results in different values of and . In general, those sensitivity coefficients generated by different methods have shown comparable results.
We assumed that a minor deviation in each input parameter, in terms of its significant digits approximation, could contribute some effects on the simulation outcome that could possibly explain the discrepancy between the simulated results and Philip’s semianalytical solution (Figure 1). In other words, the parameter values in terms of significant digits approximation that were used in computer simulation by Haverkamp et al. [7] could be different from the exact data, in terms of input parameter significant digits, that they published. Thus, we take advantage on the positive and negative relations generated from the sensitivity analysis and set up a hypothetical approximation value in Table 3 for further investigation. The cumulative effect was studied by manipulating an input parameter used for each simulation, and the subsequent manipulation of input parameter was carried out on top of the previous changed input parameter. This process begins from step 1 for base case to step 10 for spatial spacing size. For instance, the (initial cond) value (0.2376499 m3 m−3) was used as a second simulation (in step 2) after the base case simulation. This was followed by third simulation (in step 3) using value as 0.124499 m3 m−3 by remaining (initial cond) value used in the second simulation. For each simulation, Equation (11) was used to calculate the discrepancy between simulation results and Philip’s semianalytical solution (data from [7]) for absolute residual error (MA). Of all those eleven parameters in Table 3, and were the only two parameters without any limit of variation, and for this reason, we extend the variation limit by reducing the former and the latter by 98% and 90% from 500 s and 1 cm to 10 s and 0.1 cm, respectively. The and values are negative and positive relations, respectively. Decreasing and increasing the corresponding former and latter values would result in simulation failure; thus, those two parameters remained unchanged.
A consistent reduction in MA value from (initial cond) to input parameter was observed, except a slight increment at input parameter simulation; and a steep slide of MA value was observed on the input parameter simulation (Figure 3). Although the sensitivity coefficient in Figure 2 indicates that reducing value should lead to a reduction in MA value, the simulated result showed an increase in the MA value. This observation could be explained from the result of elementary effects method. This was because has different values of and , which indicate the capability of sign oscillation (Table 2).
Figure 3 shows that the simulation on the cumulative effect of steps 2–9, which combined the effect from (initial cond) (step 2) with (step 9), did not contribute to any significant effects on the advancement of water infiltration front. It only resulted in a reduction of 7.8% in MA value, from 0.0254 to 0.02343 m3 m−3. In addition, those eight input parameters had to vary in significant digits approximation as tabulated in Table 3 in order to result in the stated percentage reduction. Therefore, the significant digits approximation might not be the main cause of the problem in considering that a more significant effect on the advance of water infiltration front was shown by in the Figure 4. A further step to include in the simulation, that is, the cumulative effect of steps 2–10, which combined the effect from (initial cond) (step 2) with (step 10), there was 54.7% reduction in MA value of step 9, from 0.02343 to 0.0106 m3 m−3. This indicates that the spatial spacing size was the main cause in the advance of water infiltration front. Therefore, the simulation was repeated for the last time for the effect of spatial spacing size, alone, and in Figure 4, there was a good agreement between the simulation results and the Philip(K). This observation could be explained using (1), after rearranging it into the following form, which we termed as percentage variation in simulation results:where is the normalized input parameter value, is the normalized output parameter value, and is the normalized sensitivity coefficient (%/%). Equation (12) is simply a multiplication of the percentage change in input parameter value from the base case and the normalized sensitivity coefficient.
Using (12), the percentage variation in simulation results from input parameters of and caused an increment of 4.95% and 0.06%, respectively, despite having the highest reduction in percentage (−98%) from base case. This observation could be summarized as follows: firstly, input parameter with the highest sensitivity coefficient does not guarantee the greatest effect on the simulation result, for example, (initial cond); secondly, input parameter with the highest percentage of change also does not guarantee the greatest effect on the simulation result, for example, ; and therefore, only the highest sensitivity coefficient with the highest percentage change on input parameter (or the uncertainty) would give the most substantial effect on simulation result.
4. Conclusions
The governing equation of transient water flow in unsaturated and nonisothermal conditions was approximated numerically by finite-difference solution. It was successfully implemented into FORTRAN programming language, simulated, and verified by Philip’s semianalytical solution on water infiltration into Yolo light clay, with data from literatures.
One-at-a-time, OAT, and elementary effects, EE, methods were used in the sensitivity analysis. A common trend of sensitivity was observed across the methods in both positive and negative relations. The latter method allowed exploration of additional characteristics of input parameters at different input space, such as linearity and sign oscillation effect. The sign oscillation effect observed on input parameters explained the possibility of its deviation from those observed in OAT method at different input spaces.
A hypothetical case that was established to study the cumulative effect of input parameters on the discrepancy between simulated result and Philip’s semianalytical solution, in terms of significant digits approximation (from base case), was found to be unlikely. A large normalized sensitivity coefficient was made with initial volumetric water content and the largest percentage changes were with time-step size, but surprisingly none of them contributes to any substantial impact on simulation results, when compared to spatial spacing size. This observation led to the conclusion that the uncertainty of input parameter and normalized sensitivity coefficient of input parameters both controlled the outcome of simulation.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to acknowledge the financial support from Meiji University, Japan, and the Ministry of Education, Malaysia, and also Universiti Malaysia Terengganu, Malaysia.