International Scholarly Research Notices

Volume 2015 (2015), Article ID 824721, 7 pages

http://dx.doi.org/10.1155/2015/824721

## Sensitivity Analysis and Validation for Numerical Simulation of Water Infiltration into Unsaturated Soil

^{1}Graduate School of Agriculture, Meiji University, 1-1-1 Higashimita, Tama-ku, Kawasaki 214-8571, Japan^{2}School of Ocean Engineering, Universiti Malaysia Terengganu, 21030 Kuala Terengganu, Terengganu, Malaysia^{3}School of Agriculture, Meiji University, 1-1-1 Higashimita, Tama-ku, Kawasaki 214-8571, Japan

Received 23 July 2015; Accepted 31 August 2015

Academic Editor: Robert Černý

Copyright © 2015 Eng Giap Goh and Kosuke Noborio. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 (m^{3} 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, (m^{3} m^{−3}) and (m^{3} 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, (m^{3} m^{−3}), (m^{3} m^{−3}), and (m^{3} 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, (m^{3} m^{−3}) is residual volumetric water content, (m^{3} 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 m^{3} m^{−3}. Lower boundary was set as free-drainage to water flow. Upper boundary was set at 0.495 m^{3} 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} m^{3} 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.