High-order averaging method of hydraulic conductivity for accurate soil moisture modeling

Journal of Hydrology xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

High-order averaging method of hydraulic conductivity for accurate soil moisture modeling q Hyunuk An a, Seong Jin Noh b,⇑ a

Department of Civil Engineering, University of Suwon, San 2-2, Wau-ri, Bongdam-eup, Hwaseong-si, Gyeonggi-do 445-746, Republic of Korea Water Resources Research Division, Department of Water Resources & Environmental Research, Korea Institute of Construction Technology, 2311, Daehwa-dong, Ilsanseo-gu, Goyang-si 411-712, Republic of Korea b

a r t i c l e

i n f o

Article history: Available online xxxx Keywords: Richards’ equation Soil moisture Averaging method Hydraulic conductivity

s u m m a r y Richards’ equation (RE) is the most common mathematical expression for soil water movement in a porous medium. Despite advancements in numerical schemes and high-performance computing, the requirements of iterative computations and fine grids hinder further extension of the RE to multi-dimensional and large-scale applications. Averaging methods of hydraulic conductivity have been known to be one of the significant factors affecting the accuracy of numerical solutions of the RE, especially when coarse grids are used. In this study, we developed a high-order averaging method of hydraulic conductivity for accurate numerical modeling of the RE, which has a straightforward formula regardless of the soil conditions and produces high simulation accuracy when used on coarse grids. The developed method is based on the high-order upwind scheme, which is widely used for hyperbolic partial differential equations within a finite volume framework in order to prevent numerical oscillations near a discontinuity while preserving high-order accuracy. Numerical simulations of several one- and two-dimensional cases performed in the study indicate that the proposed method outperforms existing simple averaging methods and is also superior, or at least equivalent, to complex averaging methods over a wide range of soil textures, especially on coarse grids. In addition, the proposed method is straightforwardly extended to nonorthogonal grids by being combined with the coordinate transformation method and the extension is verified through multi-dimensional test cases as well as tests on a heterogeneous soil domain. Ó 2013 The Authors. Published by Elsevier B.V. All rights reserved.

1. Introduction Modeling of soil water movement is of great importance in various research areas, and significant advances have been made toward understanding soil moisture dynamics in the last decade (Bárdossy et al., 1995; Bronstert and Plate, 1997; Corradini et al., 2011, 1998; DiCarlo, 2013; among others). Soil moisture influences a range of environmental processes in a nonlinear manner, e.g. the partitioning of precipitation into runoff, infiltration, leakage, and evaporation; the mass and energy exchanges between soil and atmosphere; and the organization of natural ecosystems and biodiversity (Ridolfi et al., 2003; Vereecken et al., 2008; Western et al., 2002; and references therein). Richards’ equation (RE) is the most common mathematical expression for soil water movement in a porous medium but does not have a general analytical solution. Despite important advances having been made in numerical

q This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. ⇑ Corresponding author. E-mail address: [email protected] (S.J. Noh).

solutions of the RE (An et al., 2012, 2011, 2010; Celia et al., 1990; Jones and Woodward, 2000; Kurázˇ et al., 2010; Lott et al., 2012), simulation of soil moisture using the RE still requires considerable computational effort, which hinders its extension to multi-dimensional and large-scale applications. Although relatively coarse grids are required to simulate a wide spatial domain for a long time period, previous studies have shown that simple conventional methods using coarse grids suffer from drawbacks of numerical instability and oscillation, resulting in large errors (e.g. Baker, 2006). For example, if the numerical accuracy obtained at a spatial resolution of 10 cm is available with a spatial resolution of 100 cm, we can simulate a 100-fold larger domain in two-dimensional space with the same numerical computation and also expand the available temporal length in a given time. However, complex and sensitive numerical conditions originating from heterogeneous and anisotropic features of the soil matrix compound the difficulties in employing coarse grids. Therefore, accurate numerical methods on relatively coarse grids are crucial for applying numerical solutions of the RE to real-world cases. One of the significant factors affecting the accuracy of numerical solutions of the RE is the averaging method of hydraulic conductivity when coarse grids are used. Previous studies (Baker,

0022-1694/$ - see front matter Ó 2013 The Authors. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

2

H. An, S.J. Noh / Journal of Hydrology xxx (2014) xxx–xxx

1995; Belfort and Lehmann, 2005; Haverkamp and Vauclin, 1979; Szymkiewicz, 2009) have shown that simple averaging methods such as arithmetic, geometric, or upwind means may produce unacceptably large errors on coarse grids owing to the high nonlinearity of the hydraulic conductivity. Further, the arithmetic mean— one of the most popular averaging methods—often suffers from the drawback of numerical oscillations on coarse grids. In order to achieve accuracy within an acceptable level on coarse grids, the Darcian mean approach (Warrick, 1991) has been proposed. In this approach, the average hydraulic conductivity is estimated using Darcy’s law by assuming a steady state between two adjacent nodes or cells. Subsequently, several methods based on the Darcian mean approach have been proposed (Baker, 2000, 1995; Gastó et al., 2002; Szymkiewicz, 2009). Most recently, Szymkiewicz (2009) proposed the Darcian-mean-type method, in which the flow is categorized into three types—infiltration, drainage, and capillary rise—and hydraulic conductivity is computed between the integrated mean and the upwind mean according to flow types. This method generally shows a good performance for a variety of soil textures and cases. However, the integrated mean, used for computation of the Darcian mean, incurs additional computational cost compared to the simple averaging method. Since the function of hydraulic conductivity is highly nonlinear and no general form exists to integrate the hydraulic conductivity function, numerical integration is usually required as an additional step. Belfort et al. (2013) showed that only the upwind mean can strictly preserve the monotonicity of the solution. They also proposed a simple switching method, in which the upwind mean is used when the arithmetic or geometric mean violates the monotonicity of the solution. Despite this method having a simple implementation and incurring lower additional computational cost, its accuracy should be determined between the arithmetic mean and the upwind mean. In addition, if the arithmetic or geometric mean violates the monotonicity in several steps, e.g. multi-dimensional cases, this method may become identical to that with the upwind mean. The objective of this study is to develop and verify a new high-order method for averaging the hydraulic conductivity to obtain accurate numerical solutions of the RE in conjunction with coarse grids. The developed method is based on the highorder upwind scheme, which has been widely used since the 1980s for hyperbolic equations such as the compressible Euler equation (e.g. Chakravarthy and Osher, 1983; Harten et al., 1983; Roe, 1986; Sweby, 1984) or the shallow water equation (e.g. An and Yu, 2012; Louaked and Hanich, 1998; Mingham and Causon, 1998; Toro, 2000) for preventing numerical oscillations while capturing discontinuity accurately. In this study, the high-order upwind scheme is implemented for evaluating hydraulic conductivity. To the best of our knowledge, the highorder upwind scheme has not yet been applied to the evaluation of hydraulic conductivity in the numerical modeling of the RE, despite its high potential to improve the numerical accuracy on coarse grids. It is worth noting that the developed method is different from the conventional upwind mean, which has only first-order spatial accuracy. The remainder of this paper is organized as follows. Section 2 presents the numerical scheme of the RE with averaging methods of hydraulic conductivity. It also describes several conventional averaging methods and the proposed method. Section 3 presents an evaluation of the proposed method through several test cases, including one-dimensional (1D) infiltration, drainage, and capillary rise problems; and two-dimensional (2D) infiltration and drainage problems. The 2D test cases consider the implementation of the proposed method on nonorthogonal grids and a heterogeneous soil domain. Lastly, Section 4 concludes the paper.

2. Numerical scheme 2.1. Numerical discretization of RE The 1D form of the RE is written as follows (Richards, 1931):

@hðwÞ @q  ¼ 0; @t @z q ¼ KðwÞ

ð1Þ

@ ðw þ zÞ; @z

ð2Þ

where w is the pressure head; h, the volumetric soil moisture content; K, the hydraulic conductivity; t, the time; and z, the vertical dimension, which is assumed to be positive in the upward direction. As Eqs. (1) and (2) include both h and w, this form of the RE is called the ‘‘mixed’’ form. This form is generally preferred to the h-based and w-based forms for its improved mass balance in treating saturated and unsaturated flows simultaneously (Celia et al., 1990). Applying implicit Euler temporal discretization and cell-centered finite-volume spatial discretization on a uniform grid to Eq. (1) gives mþ1 mþ1 mþ1 mþ1 K mþ1 Þ K mþ1  wmþ1 K mþ1 hmþ1  hm i1 Þ iþ1=2 ðwiþ1  wi i1=2 ðwi iþ1=2  K i1=2 i i   þ Dt D z2 D z2 Dz

! ¼ 0; ð3Þ

where the superscript m denotes the time level and the subscript i denotes the cell number. Note that fourth term on the left-hand side of Eq. (3) represents nonlinear advection originating from gravity force, and this term can cause numerical oscillation near the discontinuity unless an appropriate averaging method of hydraulic conmþ1 ductivity for K mþ1 iþ1=2 and K i1=2 is implemented. Similarly, a multidimensional RE on an orthogonal grid can also be discretized within the cell-centered finite-difference or finite-volume framework, as shown in the previous studies (An et al., 2011; Clement et al., 1994). A discretized equation is solved using the Newton iteration method. Since this method is usually sensitive to the initial guess, we implement a simple line-search approach to improve the robustness following Kelley (1995). The time-step duration is adjusted on the basis of the number of iterations required for convergence at the previous time-step duration (An et al., 2012; Paniconi and Putti, 1994). The time-step duration cannot be less than a preselected minimum or greater than a preselected maximum. If the number of iterations to convergence is less than Nm, the next time-step duration is multiplied by Cm, which is a predetermined value greater than 1. If the number of iterations is greater than Nr, the next time-step duration is multiplied by Cr, which is a preselected value less than 1. Further, if the number of iterations becomes greater than a prescribed Nb, the iterative process for the time level is terminated. Subsequently, the time-step duration is multiplied by Cb—a predetermined value less than 1—and the iterative process restarts. For most large and difficult problems, these control factors of the time-step duration will require adjustment to ensure good iterative performance. In this study, the values of Cm = 1.2, Cr = 0.8, Cb = 0.5, Nm = 6, Nr = 10, and Nb = 20 are used according to the previous study (An et al., 2012), where the parameters are optimized empirically. As we deal with the RE in the ‘‘mixed’’ form, soil water retention and hydraulic conductivity models are required to express the relationships of h–w and K–w, which impose difficulties in the numerical computation owing to high nonlinearity. Several mathematical models have been proposed in the literature (Szymkiewicz and Helmig, 2011; and references therein). In this study, because of a wide range of applications, we use two sets of soil water retention and hydraulic conductivity models: (1) the combination of the Brooks–Corey model (Brooks and Corey, 1964) and Mualem model

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

3

H. An, S.J. Noh / Journal of Hydrology xxx (2014) xxx–xxx

(Mualem, 1976), and (2) the van Genuchten–Mualem model (Mualem, 1976; van Genuchten, 1980). The combination of the Brooks– Corey model (Brooks and Corey, 1964) and the Mualem model (Mualem, 1976) gives

h  hr Se ¼ ¼ hs  hr

(

1

)

ðajwjÞk

ð4Þ

;

K ¼ K s S2:5þ2=k ; e

ð5Þ

where Se is the effective saturation; hr and hs are the residual and saturated water contents, respectively; a and k are parameters whose values depend on the soil properties; and Ks is the saturated hydraulic conductivity. On the other hand, the van Genuchten– Mualem model (Mualem, 1976; van Genuchten, 1980) gives

Se ¼

h  hr ¼ hs  hr



1 1 þ ðajwjn Þ

11=n

ð6Þ

;

  11=n 2 K ¼ K s S1=2 1  1  Sn=ðn1Þ ; e e

ð7Þ

where a and n are parameters. It is noteworthy that the two selected sets of models have higher nonlinearity than other models, which also provide appropriate benchmark conditions for comparing averaging methods for hydraulic conductivity in the Section 3. 2.2. Conventional methods for averaging hydraulic conductivity As presented in previous studies (Belfort et al., 2013; Szymkiewicz, 2009), the averaging method of hydraulic conductivity at the cell interface (e.g. Ki+1/2) can significantly affect the accuracy of the numerical model on relatively coarse grids. The conventional popular averaging methods are expressed as follows:

Arithmetic mean : K arit ¼ 0:5ðK iþ1 þ K i Þ; Geometric mean : K geo ¼

ð8aÞ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi K iþ1 K i ;

ð8bÞ

Harmonic mean : K har ¼ ð1=ð1=K iþ1 þ 1=K i ÞÞ; Integrated mean : K int ¼  Upwind mean : K up ¼

1 wiþ1  wi K iþ1 ; K i;

Z

ð8cÞ

wiþ1

KðwÞdw;

ð8dÞ

wi

2.3. Proposed high-order method In this study, we develop a high-order method, in which the average hydraulic conductivity is determined by a straightforward formula regardless of the soil moisture conditions. The proposed method is described as follows. In the cell-centered finite-volume or finite-difference model, the value of K is defined at the cell center. Piecewise linear distribution of K is assumed inside the cell. As a result, two different values of K exist at the cell interface (i + 1/2), as shown in Fig. 1. Khigh is selected in an upwind manner as

( K high ¼

K iþ1=2þ ;

if /iþ1  /i

K iþ1=2 ;

else

ð9Þ

;

where Ki+1/2+ and Ki+1/2 are reconstructed values on the upper and lower sides, respectively, of the cell interface (i + 1/2), and are calculated by

K iþ1=2þ ¼ K iþ1  ðrK iþ1  Dziþ1 Þ=2;

ð10Þ

K iþ1=2 ¼ K i þ ðrK i  Dzi Þ=2;

where rKi is the gradient of K in cell i. Note that if the gradient of K is zero, this method becomes identical to the conventional upwind mean, which results in first-order spatial accuracy. In this study, the gradient of K is evaluated as

rK i ¼ /ðri ÞDK i1=2 ;

if /iþ1  /i ; else

ð11Þ

ð8eÞ

where /i = wi + zi is the total head. Previous studies have shown that none of the above-listed averaging schemes produce satisfactory results for a wide range of soil textures and flow directions on coarse grids (Szymkiewicz, 2009). With the aim of achieving generally accurate results on coarse grids, several Darcian mean approaches have been proposed (Baker, 2000, 1995; Gastó et al., 2002; Szymkiewicz, 2009). Among those, the most recently proposed Darcian mean is given as (Szymkiewicz, 2009)

K szy

Note that the direction of z is opposite in Szymkiewicz (2009). This method shows a generally satisfactory performance for a wide range of soil textures and situations. Although this algorithm is rather complex, the additional cost incurred is not significant, except for the cost for computing Kint. To avoid incurrence of a large computational cost for the numerical integration of K, Szymkiewicz (2009) generated a precomputed table for Kint before the simulation, and then interpolated the value of Kint using this table. This approach is efficient for a homogeneous or layered soil domain. However, its application to a highly heterogeneous soil domain is unfeasible owing to its memory constraint, because the table for Kint should be generated for each type of soil texture. It is beyond the scope of this study to provide a computationally efficient way for the Darcian mean. Therefore, Kint is just computed at every iteration step by numerical integration.

  8 K > max K int ; 1Diþ1 ; if Dw=Dz < 0 > w=Dz > > <   K iþ1 2 ¼ min 1Dw=Dz ; Kðwi  Dw =DzÞ ; else if 0  Dw=Dz < 1 > > > > DzK ð1Þ K ð2Þ : ; else ðDzdzÞK ð1Þ þdzK ð2Þ

z

i+1

K i+1 K i+1/2+

i

K i+1/2Ki

i-1

K i-1

ð8fÞ where Dw = wi  wi+1, K(2) = K(wi  Dz) and

dz ¼

 Dw þ

Dz = zi+1  zi,

K(1) = Kint(wi+1, wi  Dz),

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dw2 þ 4ðK ð1Þ =K ð2Þ  1ÞðDw  DzÞDz 2ðK ð1Þ =K ð2Þ  1Þ

K :

Fig. 1. A schematic of second-order upwind scheme. Piecewise linear distribution of hydraulic conductivity is assumed.

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

4

H. An, S.J. Noh / Journal of Hydrology xxx (2014) xxx–xxx

where ri = DKi+1/2/DKi1/2, DKi1/2 = Ki  Ki1, DKi+1/2 = Ki+1  Ki, and /(r) is a slope limiter function, which prevents local extrema at the cell interface, thus preserving monotonicity. Several choices are available for the slope limiter (Roe, 1986; Sweby, 1984); the minmod limiter is employed in this study for its superior numerical stability. This limiter is expressed as

/ðrÞ ¼ maxð0; minðr; 1ÞÞ:

ð12Þ

Note that the values of Ki1, Ki, and Ki+1 are utilized for gradient evaluation in the cell i. The proposed method is a modified version of the high-order upwind scheme. Although the RE is not a hyperbolic equation, it has an advection-like term originating from the gravitational force, as well as second-order derivative terms. However, the high-order upwind scheme is usually used in combination with explicit temporal discretization. In this study, for mitigating the convergence problem in the highly nonlinear RE, which requires implicit temporal discretization and Newton-like iteration, the high-order upwind scheme is implemented only for computing the hydraulic conductivity at the cell interface. 2.4. Numerical treatment for boundary condition The averaging method of hydraulic conductivity can significantly affect the numerical accuracy under the Dirichlet boundary condition, also known as the fixed boundary condition. The implementation of the proposed method under the Dirichlet boundary condition is not straightforward, because insufficient information is available at the boundary interface. As a result, the accuracy at the Dirichlet boundary can deteriorate. Then, a different treatment is proposed in this study for the Dirichlet boundary condition. If the Dirichlet boundary value of wb, where subscript b denotes the boundary, is set at the boundary interface in cell i, the flux through the boundary interface is given as

q ¼ K av

  @w þ1 ; @z

ð13Þ

where Dwb = wb  wi, Dzb = zb  zi, zb is the z value at the boundary, and Kav is the averaged hydraulic conductivity between the cell center and the boundary interface. Note that q in Eq. (13) is added on the right-hand side of Eq. (3) as a source term. Here, we estimate the flux through the boundary interface as

! @w þ1 ; @z up

q ¼ K up

ð14Þ

where the gradient of the pressure head is assumed to be estimated on the upwind side as

@w @z

¼

up

@w @K ; @K up @z up

ð15Þ

where

@K Kðwb Þ  Kðwi Þ ¼ ; @z up Dz b

ð16Þ

and

@w dw ¼ ; @K up Kðwup Þ  Kðwup  dwÞ

ð17Þ

where wup is the pressure head on the upwind side, and is expressed as

 wup ¼

wb ;

if wb þ zb  wi þ zi

wi ;

else

;

and dw is defined according to a given case as

ð18Þ

8 if  Dwb =Dzb < 0 > < eDwb ; dw ¼ Dwb ; else if 0 < Dwb =Dzb  1 ; > : minðDwb ; Dzb Þ: else

ð19Þ

where e is a small dimensionless number, which is set as e = 103 in this study. The cases considered in Eq. (19) are infiltration, drainage, and capillary rise, respectively, which are the same as those in Eq. (8f).

3. Numerical evaluation The performance of the proposed method is verified by performing several numerical evaluations. The simulation conditions of test cases are listed in Table 1, including the dimensions, soil moisture conditions, soil models, texture (Carsel and Parrish, 1988; Rawls et al., 1991), simulation periods, initial conditions, boundary conditions, spatial resolutions, and flow domains. The test cases are designed to verify the performance of the proposed method for various soil moisture conditions as well as for a relatively wide range of soil textures, whose property parameters are listed in Table 2. Therefore, each test case has different setups for simulating relevant soil moisture conditions effectively. For example, the flow domain for capillary rise is 0.5 m because capillary rise will be simulated only near the surface, whereas other 1D cases have a flow domain of 5 m. A set of soil–water retention and hydraulic conductivity models is implemented in turn on a homogeneous soil domain for three soil textures (1D) and on heterogeneous soil domain (2D). It is also worthy to note that three soil textures in one dimension are selected to cover various soil types and that the performance results with the other soil textures having relatively similar saturated hydraulic conductivity (e.g. sandy loam and loam, or silty loam and clay) are not significantly different (not shown). 3.1. One-dimensional cases The 1D test cases include infiltration, drainage, and capillary rise problems; the simulation of each test case is performed on a homogeneous soil domain with three typical soil textures in order to verify the performance of the proposed method for a wide range of soil textures. For comparison purposes, four existing averaging methods, i.e., Kari, Kup, Kint, and Kszy, are also used in these test simulations. Khar and Kgeo are not considered in this study, because their insufficient behaviors have been well reported in other papers (Belfort et al., 2013; Szymkiewicz, 2009) and treating too many methods hinders a clear description of the simulation results. For quantitative analysis of error, two error measurement metrics are evaluated with mesh refinement following Hu et al. (2012) and Szymkiewicz and Helmig (2011):

Eq ¼

jqRef  qi j ; jqRef j

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2ffi u N u w  w 1X i t Ref ; Ew ¼ N i¼1 wRef

ð20Þ

ð21Þ

where qRef and qi are the volumes of infiltrated, exfiltrated, and evaporated water at the boundary obtained from the reference solution and from the solution of the different averaging methods, respectively; wRef and wi are the reference pressure head and the pressure head obtained from the different averaging methods, respectively; and N is the number of cells.

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

5

H. An, S.J. Noh / Journal of Hydrology xxx (2014) xxx–xxx Table 1 Simulation conditions of test cases. 1D/2D

Condition Soil model Soil textures Simulation period

Initial condition Boundary condition

Spatial resolution

Flow domain

1D

Infiltration Van Sand Genuchten Sandy loam Silty loam Drainage Brooks– Sand Corey Loam Clay Capillary Van Sand rise Genuchten Loam Clay

winit = 5 m

Dz = 5 10, 25, 50, 100 cm

5 m homogeneous soil domain

Dz = 5 10, 25, 50, 100 cm

5 m homogeneous soil domain

Dz = 0.5, 1, 2.5, 5, 10 cm

0.5 m homogeneous soil domain

Dz = 5 10, 25, 50, 100 cm

5 m  5 m heterogeneous soil domain 5 m  5 m layered soil domain

2D

2h 14 h 120 h 30 h

wb = 0.01 m at z=0m

winit = 0.075 m wb = 0.075 m at z = 5 m

7d

winit = 0.5–z m wb = 10 m at z=0m

Infiltration Van Sand + Loam 120 h Genuchten Drainage Brooks– Sand + Loam 30 h Corey

winit = 5 m

wb = 0.01 m at

z=0m winit = 0.075 m wb = 0.075 m at z = 5 m

Non-uniform, 5  5, 10  10, 20  20, 50  50

Table 2 Soil property parameters used in test simulations. Model

Soil texture

hS (m3/m3)

hr (m3/m3)

a (m1)

n or k

Reference

5

KS (m/s)

Van Genuchten model

Sand Sandy loam Loam Silty loam Clay

0.43 0.41 0.43 0.45 0.38

0.045 0.065 0.078 0.067 0.068

8.250  10 1.228  105 2.889  106 1.250  106 5.556  107

14.5 7.5 3.6 2.0 0.8

2.68 1.89 1.56 1.41 1.09

Carsel and Parrish (1988)

Brooks–Corey model

Sand Loam Clay

0.44 0.46 0.38

0.020 0.027 0.068

5.830  106 3.667  106 1.667  108

13.8 9.0 2.7

0.59 0.22 0.13

Rawls et al. (1991)

3.1.1. Infiltration An infiltration problem is imposed on a 1D vertical domain (depth: 5 m). Three types of soil textures are considered: sand, sandy loam, and silty loam of the van Genuchten model listed in Table 1. The initial condition is set as winit = 5 m, and the Dirichlet boundary of wb = 0.01 m is set at the top of the flow domain, which imposes the dry initial condition and an almost saturated boundary condition, making the numerical computation more challenging on account of the high nonlinearity of hydraulic conductivity. The simulations are performed on five different grids, as Dz = 5, 10, 25, 50, 100 cm. The simulation periods for sand, sandy loam, and silty loam soil textures are 2, 14, and 120 h, respectively, given that the speeds of infiltration differ according to the soil textures. Fig. 2(a) shows the pressure head profiles with a spatial resolution of Dz = 50 cm at the end of simulations of the three textures. Since the integrated mean has a convergence problem in the simulation of silty loam, the result of Kint is not shown in the silty loam simulation result in Fig. 2(a). ‘‘Ref’’ indicates the reference solution, which is obtained by Kup on a dense grid. Before comparison of the results, it was confirmed that most averaging methods gave very similar results on a dense grid. For all soil textures except silty loam, the position of the wetting front obtained by the proposed method is quite similar to that obtained by Szymkiewicz’s method. In the case of silty loam, the proposed method captures the position of the wetting front better than does Szymkiewicz’s method. The arithmetic mean generates numerical oscillation in the simulations of sand and sandy loam, but this is not observed for silty loam. Overall, the integrated mean underestimates the infiltration and the upwind mean overestimates it. These results are consistent with the simulation results of Belfort et al. (2013) and Szymkiewicz (2009). The estimated errors in the simulation of infiltration with different averaging methods and spatial resolutions are shown in Fig. 2(b). It is clear that the proposed method shows the best

performance for all three textures. After the proposed method, the arithmetic mean produces superior results to the upwind mean and Szymkiewicz’s method on fine grids, but it often produces larger errors on coarse grids. In contrast, Szymkiewicz’s method is generally superior to the other conventional averaging methods on relatively coarse grids but the errors of this method converge slowly because this algorithm converges to the integrated mean in all cases of infiltration when Dz is close to zero. The integrated mean seems to be inappropriate for the simulation of infiltration. Note that Szymkiewicz’s method shows a comparatively larger error in the simulation of silty loam. For investigating the errors of different averaging methods on coarse grids in greater detail, the volumes of infiltrated water with Dz = 50 cm over time were obtained from different averaging methods except for the integrated mean, given that the integrated mean produces poor results as shown in Fig. 2(b); the results are shown in Fig. 2(c). It is interesting to note that the arithmetic mean underestimates or overestimates the infiltration depending on the soil texture, whereas the upwind scheme always overestimates it. In addition, Szymkiewicz’s method shows a constant infiltration rate when Dz = 50 cm. On relatively coarse grids, Szymkiewicz’s method seems to maintain a steady infiltration rate at the boundary. This method is effective for soils having a highly nonlinear function of hydraulic conductivity, such as sand or sandy loam; however, it fails to simulate the infiltration on soils having a mildly nonlinear function of hydraulic conductivity, such as silty loam. In contrast, the proposed method approaches the actual infiltration volume with time for all soil textures. 3.1.2. Drainage The second test case is the 1D drainage problem following Szymkiewicz (2009). In this case, the domain is considered to be 5 m in depth. The initial condition is set as winit = 0.075 m, and the lower boundary is maintained at the initial pressure head. The simulation period is 30 h. Three types of soil textures are

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

H. An, S.J. Noh / Journal of Hydrology xxx (2014) xxx–xxx

0

0

-0.5

-0.5

-1

-1

-1.5

-1.5

-2

z (m)

z (m)

(a)

-2.5 -3

0

-1

-2

z (m)

6

-2.5 -3

Ref Kari Kup Kint Kszy Khigh

-3.5 -4

sand

-4.5 -5 -5

-4

-3

-2

sandy loam

-4 -4.5 -5

-1

0

1

-5

-4

-3

Kari Kup Kint Kszy Khigh

Ref Kari Kup Kszy Khigh

silty loam

-4

-5

-1

0

1

-5

-4

-3

-1

0

Kari Kup Kszy Khigh 10

10

Eq (%)

Eq (%) 1

-2

ψ (m)

Kari Kup Kint Kszy Khigh

100

Eq (%)

10

-2

-3

ψ (m)

ψ (m)

(b) 100

Ref Kari Kup Kint Kszy Khigh

-3.5

-2

1 1

sandy loam

sand 0.1

0.1

10

10

100

10

100

Δz (cm)

0.7

sand

0.5 0.4 0.3 0.2

Ref Kari Kup Kszy Khigh

0.1 0

0

0.5

1

1.5

2

sandy loam

0.6

Infiltrated water volume (m)

0.6

Infiltrated water volume (m)

Infiltrated water volume (m)

100

Δz (cm)

Δz (cm)

(c)

silty loam

0.5 0.4 0.3 Ref Kari Kup Kszy Khigh

0.2 0.1 0

0

2

time (hr)

4

6

8

time (hr)

10 12 14

0.5

silty loam

0.4 0.3 0.2 Ref Kari Kup Kszy Khigh

0.1 0

0

20

40

60

80 100 120

time (hr)

Fig. 2. (a) The pressure head profiles obtained from different averaging methods with Dz = 50 cm at the end of simulation, where Ref is the reference solution, (b) the errors of infiltration obtained from different averaging methods, and (c) infiltrated water volumes obtained from different averaging methods with Dz = 50 cm in 1D infiltration test cases.

considered in this case: sand, loam, and clay of the Brooks–Corey model listed in Table 2. The simulations are performed on five different grids, with Dz = 5, 10, 25, 50, 100 cm. Fig. 3(a) shows the pressure head profiles with Dz = 50 cm at the end of simulations for the three textures. As was reported in previous studies (Belfort et al., 2013; Szymkiewicz, 2009), in this study, too, the arithmetic and integrated means caused numerical oscillations in the simulations of clay as well as sand. For all three textures, no numerical oscillation is observed in the simulation using the other three methods. However, the accuracies of solutions obtained from the upwind mean and Szymkiewicz’s method are degraded near the surface of the domain (z > 1 m), but the solutions obtained from the proposed method are quite close to the reference solutions near the surface of the domain in all cases. The estimated errors in the simulation of exfiltration with the different averaging methods are shown in Fig. 3(b). The proposed method shows the best performance for sand and loam; the integrated mean shows good performance for clay but causes numerical oscillation in the case of sand. The proposed method may be second best for clay. Overall, the convergence speeds of errors for the five different methods seem to be similar, except for Szymkiewicz’s method. Szymkiewicz’s method shows a smaller error than

the arithmetic mean on coarse grids but shows a larger error than the arithmetic mean on relatively dense grids. The estimated errors in the pressure head with the different averaging methods are shown in Fig. 3(c). As was the case with the errors of exfiltration, the convergence speeds of errors in the pressure head for the five different methods seem to be similar. However, the proposed method produces the smallest error on coarse grids; nevertheless, on relatively dense grids, it does not show a smaller error than the other methods. As mentioned in Section 1, our key focus in this study is on the accuracy when a relatively coarse grid is used, because the accuracy of the numerical model of the RE is less affected by the averaging method of hydraulic conductivity on a fine grid. It is also found that the error behaviors of the averaging methods for drainage cases are different from those for the infiltration cases. The reason is supposedly the different nonlinearity of K depending on the case. 3.1.3. Capillary rise The third test case is the 1D capillary rise problem, which mostly corresponds to Case D in the study of Szymkiewicz and Helmig (2011); the only difference is that in our study, this test

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

7

H. An, S.J. Noh / Journal of Hydrology xxx (2014) xxx–xxx 0

0

sand

-0.5

-1

-1

-1.5

-1.5

-1.5

z (m)

-2

-2.5

-3.5 -4 -4.5

-2

-2.5

-3

-3

Ref Kari Kup Kint Kszy Khigh

-3.5 -4 -4.5

-5 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

-3.5 -4 -4.5

-5 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1

10

1

Kari Kup Kint Kszy Khigh

10

1

0.1

100

10

Kari Kup Kint Kszy Khigh

0.1

0.01

0.01

100

Δz (cm)

Kari Kup Kint Kszy Khigh

0.1

Kari Kup Kint Kszy Khigh

0.01

E ψ (m)

0.1

10

100 Δz (cm)

E ψ (m)

E ψ (m)

clay 0.01

Δz (cm)

(c)

1

loam 0.01

10

Kari Kup Kint Kszy Khigh

0.1

sand 0.1

0

ψ (m)

Eq (%)

Kari Kup Kint Kszy Khigh Eq (%)

Eq (%)

10

Ref Kari Kup Kint Kszy Khigh

ψ (m)

ψ (m)

(b)

-2.5 -3

Ref Kari Kup Kint Kszy Khigh

-5 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0

0

clay

-0.5

-1

-2

z (m)

0

loam

-0.5

z (m)

(a)

0.001 0.001

0.001

sand 0.0001

10

0.0001

100

Δz (cm)

loam 10

100

Δz (cm)

clay 0.0001

10

100

Δz (cm)

Fig. 3. (a) The pressure head profiles obtained from different averaging methods with Dz = 50 cm at the end of simulation, where Ref is the reference solution, (b) the errors of exfiltration, and (c) the errors of pressure head obtained from different averaging methods in 1D drainage test cases.

performed is on a homogeneous soil domain whereas in their study, the test is performed on a two-layered soil domain. In this study, a 0.5-m domain is considered. The initial condition is set as winit = 0.5–z m; the Dirichlet boundary of wb = 10 m is set at the top of the flow domain; and the lower boundary is set to be a no-flow boundary. Since the dry condition is maintained at the top of the flow domain, the soil water moves upward and evaporation occurs at the top of the flow domain. The capillary rise problem requires a relatively longer time to occur than the infiltration or drainage problem, and the simulation period is 7 days. Three types of soil textures are considered in this case: sand, loam, and clay of the van Genuchten model listed in Table 2. The simulations are performed on five different grids, with Dz = 0.5, 1, 2.5, 5, 10 cm. Fig. 4(a) shows the pressure head profiles with Dz = 5 cm at the end of simulation for the three textures. Since Szymkiewicz’s

method has a convergence problem in the case of clay, the result of Kszy is not shown here. The integrated mean shows the best performance for all cases, and Szymkiewicz’s method shows results that are very similar to those of the integrated mean. For sand, the proposed method produces profiles to similar to those for the integrated mean, whereas for the other two textures, the proposed method shows less accurate results than the integrated mean. However, the results obtained from the proposed method are still much better than those obtained from the upwind and arithmetic means. Overall, the upwind and arithmetic means overestimate the volume of evaporated water. The estimated errors in the simulation of evaporation with the different averaging methods are shown in Fig. 4(b). The proposed method shows the best performance for sand, whereas the integrated mean shows the best performance for the other two textures. In fact, the integrated mean is effective for the capillary

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

8

H. An, S.J. Noh / Journal of Hydrology xxx (2014) xxx–xxx

(a)

0

0

0

sand

loam

-0.1

-0.1

-0.2

-0.2

clay

-0.05 -0.1

-0.3 Ref Kari Kup Kint Kszy Khigh

-0.4

-0.5 -1

z (m)

z (m)

z (m)

-0.15

-0.3

-0.3 Ref Kari Kup Kint Kszy Khigh

-0.4

-0.5

-0.8 -0.6 -0.4 -0.2

0

-2

-0.35 -0.4 -0.45

-1.5

ψ (m)

-0.5

-0.5 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1

0

sand

100

100

ψ (m) 100

Kari Kup Kint Kszy Khigh

loam 10

Eq (%)

Kari Kup Kint Kszy Khigh

-1

Ref Kari Kup Kint Khigh

ψ (m)

Eq (%)

Eq (%)

(b) 1000

-0.2 -0.25

10

Kari Kup Kint Khigh

clay

1

0.1 1

10

0.01 1

10

1

Δz (cm)

10

1

Δz (cm)

10

Δz (cm)

Fig. 4. (a) The pressure head profiles obtained from different averaging methods with Dz = 5 cm at the end of simulation, where Ref is the reference solution, and the errors of evaporation obtained from different averaging methods in 1D capillary rise test cases.

rise simulation, as reported by Szymkiewicz (2009). Szymkiewicz’s method shows errors similar to those of the integrated mean but it causes a convergence problem in the case of clay. Although the proposed method does not produce smaller errors than those for the integrated mean in simulations of any of the textures, it still shows much better accuracy than the simple conventional averaging methods, i.e., the upwind and arithmetic means. The upwind and arithmetic means produce extremely large errors for sand on coarse grids, which is consistent with the results reported in previous studies (Szymkiewicz and Helmig, 2011; Szymkiewicz, 2009). It can be noticed that the errors for the capillary rise cases are different from those for the infiltration and drainage cases. In these cases, the soil water moves upward and the diffusion terms are dominant over the advection term originating from gravity force. Thus, these results can be translated as follows: the integrated mean is accurate when the diffusion terms are dominant, which is consistent with the observation that the integrated mean performs well for horizontal flow (Szymkiewicz, 2009), where no advection term is considered. 3.2. Multi-dimensional cases The proposed method can be extended to multi-dimensional modeling in a straightforward manner. The method can also be applied to nonorthogonal grids by being combined with the coordinate transformation method. The coordinate transformation method for the RE is briefly introduced here for completeness. A curvilinear coordinate system in physical space, (x, y, z) = (x1, x2, x3), where x and y denote horizontal dimensions, can be transformed into a new coordinate system in uniform computational space, (n, g, f) = (n1, n2, n3), by the coordinate transformation method. The summation form of 3D RE without anisotropy is written as

  3 @h X @ @ðw þ zÞ ¼ 0;  K @t r¼1 @xr @xr

ð22Þ

and the coordinate-transformed RE is given as (An et al., 2010)

  3 X 3 1 @h X @ @ðw þ zÞ ¼ 0;  Gp;q K J @t p¼1 q¼1 @np @nq

ð23Þ

where

J¼

@ðn1 ; n2 ; n3 Þ ; @ðx1 ; x2 ; x3 Þ

Gp;q ¼

3 X K @np @nq : J @xr @xr r¼1

ð24Þ

Here, J is the Jacobian determinant that represents the ratio of the control volume in physical space to that in computational space and Gp,q is an operator that represents the mesh sleekness. Since the coordinate-transformed equation, Eq. (23), is on the grid defined in computational space, in which the grid is set to be uniform, the numerical scheme for the orthogonal grid can be applied without any modification. The hydraulic conductivity in the control volume can be assumed to have a piecewise linear distribution for each of the directions of n, g, or f, and the proposed method can be implemented simply. More details of the coordinate-transformed model of the RE can be found in An et al. (2010). The performance of the proposed method is investigated for a heterogeneous soil domain in multi-dimensional test cases. The integrated mean and Szymkiewicz’s method are not considered in this investigation, because their implementations on nonorthogonal grids and heterogeneous soil are not straightforward. We note that Szymkiewicz’s method has been extended to heterogeneous soil by Szymkiewicz and Helmig (2011), but doing so in our multi-dimensional test simulations increased the nonlinearity of the system and often caused divergence during the Newton iteration. Since a thorough investigation of the Darcian mean is beyond the scope of this study, it is not presented in this paper. 3.2.1. Two-dimensional infiltration on heterogeneous soil domain Infiltration on a 2D heterogeneous soil domain is considered in this test simulation. A sand soil block (2.5 m  2 m) is located inside a square loam soil domain (5 m  5 m), as shown in Fig. 5. The van Genuchten model is used in this test; the soil parameters are listed in Table 2. The simulation settings are similar to those

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

9

H. An, S.J. Noh / Journal of Hydrology xxx (2014) xxx–xxx 0

0

K ari, Δz=50cm

Loam

-1

Sand

Sand

Z (m)

θ=0.1 θ=0.2 θ=0.3

-2

Z (m)

θ=0.1 θ=0.2 θ=0.3

-2

-3

-3

-4

-4

-5

K up, Δz=50cm

Loam

-1

-5 0

1

2

3

4

5

0

1

2

X (m) 0

-1

5

Loam

Reference solution

Sand

θ=0.1 θ=0.2 θ=0.3

-1

Sand

θ=0.1 θ=0.2 θ=0.3

-2

-2

Z (m)

Z (m)

4

0

K high, Δz=50cm

Loam

-3

-3

-4

-4

-5

3

X (m)

0

1

2

3

4

-5

5

0

1

2

X (m)

3

4

5

X (m)

Fig. 5. Contour lines of water content obtained from different averaging methods with Dx = Dz = 50 cm at the end of simulation, where Ref is the reference solution, in 2D infiltration on heterogeneous soil domain test.

(a)

(b)

0

-0.5 -1

Ref K ari K up K high

-1.5

-2

z (m)

-2

-2.5

-2.5

-3

-3

Δz=50cm

-3.5

-4

-4.5

-4.5 0.15

0.2

0.25

0.3

Δz=25cm

-3.5

-4

-5 0.1

Ref K ari K up K high

-1

-1.5

z (m)

0

-0.5

0.35

0.4

-5 0.1

0.45

0.15

0.2

Water content

(c)

0

-0.5 -1

0.25

0.3

0.35

0.4

0.45

Water content

(d)

Ref K ari Kup Khigh

10

K ari K up K high

-1.5

E q (%)

z (m)

-2 -2.5 -3

1

Δz=10cm

-3.5 -4

0.1

-4.5 -5 0.1

0.15

0.2

0.25

0.3

0.35

Water content

0.4

0.45

1

10

Δz (cm)

Fig. 6. The water content profiles obtained from different averaging methods with (a) Dx = Dz = 50 cm, (b) Dx = Dz = 25 cm and (c) Dx = Dz = 10 cm at the end of simulation, where Ref is the reference solution, and (d) the errors of infiltration obtained from different averaging methods in infiltration 2D on heterogeneous soil domain test case.

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

10

H. An, S.J. Noh / Journal of Hydrology xxx (2014) xxx–xxx 0

0

5X5 grid

10X10 grid -1

-2

-2

z (m)

z (m)

-1

-3

-3

-4

-4

-5

-5 0

1

2

3

4

5

0

1

2

x (m)

3

4

5

3

4

5

x (m) 0

0

20X20 grid

50X50 grid -1

-2

-2

z (m)

z (m)

-1

-3

-3

-4

-4

-5

0

1

2

3

4

5

-5

0

1

x (m)

2

x (m)

Fig. 7. Crumpled grids used in 2D drainage on non-orthogonal grids test cases.

(a)

0

(b) 10

-1

K ari K up K high

Eq (%)

z (m)

-2

Loam

1

-3

Sand -4

-5 -0.6

Ref K ari K up K high -0.5

0.1 -0.4

-0.3

-0.2

-0.1

ψ (m)

10

100

Δz (cm)

Fig. 8. (a) the pressure head profiles obtained from different averaging methods with 10  10 grid at the end of simulation, where Ref is the reference solution, and (b) the errors of exfiltration obtained from different averaging methods in 2D drainage on non-orthogonal grids test case.

outlined in Section 3.1.1. The initial condition is set as winit = 5 m, and the Dirichlet boundary of wb = 0.01 m is set at the top of the flow domain. The simulation period is 120 h. The simulations are performed on five different grids, with Dx = Dz = 5, 10, 25, 50, 100 cm. Fig. 5 shows the contour lines of volumetric water content obtained by the different averaging methods with Dx = Dz = 50 cm at the end of simulations. The positions of the wetting front as obtained by the upwind and arithmetic means are lower than the reference positions, but the positions obtained by the proposed method are similar to the references positions. In addition, the arithmetic mean shows numerical oscillation in the sand soil block. This oscillation can be seen more clearly in Fig. 6(a)–(c), which

illustrate the vertical water-volume profiles at x = 2.5 m as obtained by the different averaging methods with Dx = Dz = 50 cm, Dx = Dz = 25 cm, and Dx = Dz = 10 cm at the end of simulations. The arithmetic mean produces numerical oscillation even on relatively dense grids near the sand–loam interface owing to underestimation of the hydraulic conductivity. In contrast, the proposed method produces the best wetting-front position among the three different averaging methods without any numerical oscillation. Fig. 6(d) shows the estimated error in the simulation of infiltration by the different averaging methods. A similar tendency to that in the case of 1D infiltration in homogeneous soil test cases is observed. The proposed method outperformed the simple averaging methods in the 2D case.

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

H. An, S.J. Noh / Journal of Hydrology xxx (2014) xxx–xxx

3.2.2. Two-dimensional drainage on nonorthogonal grids Next, the drainage problem on a 5-m square soil domain is evaluated with crumpled grids to demonstrate the applicability of the proposed method to nonorthogonal grids. Crumpled grids, shown in Fig. 7, are generated by randomly moving the vertex nodes by Dz/4. A heterogeneous soil layer representing the sand texture of soil is set below z = 3 m and the loam soil texture is set elsewhere. The Brooks–Corey model is used, and the soil parameters are as listed in Table 2. The simulation settings are similar to those outlined in Section 3.1.2. Just like in the 1D drainage case, the initial condition is set as winit = 0.075 m, and the lower boundary is maintained at the initial pressure head. The simulation period is 30 h. The simulations are performed on four different grids—i.e., 5  5, 10  10, 20  20, 50  50 grids—as shown in Fig. 7. Fig. 8(a) shows the pressure head profiles with the 10  10 grid at the end of simulation. As was observed for the simulations of the 1D cases, the arithmetic mean produced numerical oscillation, and the upwind mean showed a less accurate pressure head profile near the surface of the flow domain. In contrast, the results of the proposed method were accurate on nonorthogonal grids with heterogeneous soil layers. The estimated errors are shown in Fig. 8(b). The errors of the proposed method converged well with mesh refinement and are smaller than those of the simple averaging methods, similar to the result for 1D drainage in the homogeneous soil case.

11

wide range of soil textures on coarse grids. Another advantage of this method is that its implementation to nonorthogonal grids in combination with the coordinate transformation method is straightforward. The need for high computational efforts with fine grids has limited the applications of the RE in a restricted spatial and temporal domain. Averaging methods of hydraulic conductivity have been known to be one of the significant factors affecting the accuracy of numerical solutions of the RE, especially when coarse grids are used. Therefore, the proposed averaging method of hydraulic conductivity is expected to contribute to extending the use of the RE to multi-dimensional and large-scale applications with coarse grids. Future study will be focused on uncertainty assessment of multidimensional soil moisture modeling in conjunction with data assimilation and accurate numerical solutions of the RE proposed in this study. Acknowledgement This work was partly supported by a grant from a Strategic Research Project funded by the Korea Institute of Construction Technology and by the Supercpmputing Center/Korea Institute of Science and Technology Information with supercomputing resources including technical support (KSC-2013-C3-031). References

4. Conclusions In this paper, we proposed a high-order averaging method of hydraulic conductivity for numerical modeling of the RE. The proposed method, modified on the basis of the high-order upwind scheme, has a straightforward formula regardless of soil conditions and produces high simulation accuracy in conjunction with coarse grids. Unlike the conventional high-order approach for a hyperbolic partial differential equation, the high-order upwind scheme is used only for evaluating the hydraulic conductivity at the cell interface. Several 1D and 2D cases were simulated numerically to verify the performance of the proposed method on homogeneous and heterogeneous soil domains with a wide range of soil textures. As 1D cases, infiltration, drainage, and capillary rise were simulated in three different homogeneous soil domains. The proposed method outperformed the conventional averaging methods in the infiltration and drainage cases, particularly on coarse grids, whereas the conventional methods produced either numerical oscillations or inaccurate results. In the capillary rise case, the integrated mean averaging method often performed better than the proposed method depending on the soil texture but the proposed method still showed much smaller error than the arithmetic and upwind means. Two-dimensional numerical simulations were performed for the infiltration case on a heterogeneous soil domain and the drainage case on nonorthogonal grids. Without further numerical treatments, the proposed method clearly outperformed the conventional simple methods on nonorthogonal grids and a heterogeneous soil domain. As stated above and also presented in Szymkiewicz (2009), the accuracy of numerical solutions for the RE is highly dependent on the averaging method of hydraulic conductivity, which may have been a motivation of several recently proposed complex or switching algorithms. However, complex forms require additional computational efforts, and several averaging methods are not straightforward in multi-dimensional applications. Although the proposed method did not always show better performance than the other conventional averaging methods, several numerical experiments showed that the proposed method works well for a

An, H., Ichikawa, Y., Tachikawa, Y., Shiiba, M., 2010. Three-dimensional finite difference saturated-unsaturated flow modeling with nonorthogonal grids using a coordinate transformation method. Water Resour. Res. 46, W11521. An, H., Ichikawa, Y., Tachikawa, Y., Shiiba, M., 2011. A new iterative alternating direction implicit (IADI) algorithm for multi-dimensional saturated– unsaturated flow. J. Hydrol. 408, 127–139. An, H., Ichikawa, Y., Tachikawa, Y., Shiiba, M., 2012. Comparison between iteration schemes for three-dimensional coordinate-transformed saturated–unsaturated flow model. J. Hydrol. 470–471, 212–226. An, H., Yu, S., 2012. Well-balanced shallow water flow simulation on quadtree cut cell grids. Adv. Water Resour. 39, 60–70. Baker, D., 2006. General validity of conductivity means in unsaturated flow models. J. Hydrol. Eng. 11, 526–538. Baker, D.L., 1995. Darcian weighted interblock conductivity means for vertical unsaturated flow. Ground Water 33, 385–390. Baker, D.L., 2000. A Darcian integral approximation to interblock hydraulic conductivity means in vertical infiltration. Comput. Geosci. 26, 581–590. Bárdossy, A., Bronstert, A., Merz, B., 1995. 1-, 2- and 3-dimensional modeling of water movement in the unsaturated soil matrix using a fuzzy approach. Adv. Water Resour. 18, 237–251. Belfort, B., Lehmann, F., 2005. Comparison of equivalent conductivities for numerical simulation of one-dimensional unsaturated flow. Vadose Zone J. 4, 1191–1200. Belfort, B., Younes, A., Fahs, M., Lehmann, F., 2013. On equivalent hydraulic conductivity for oscillation–free solutions of Richard’s equation. J. Hydrol. 505, 202–217. Bronstert, A., Plate, E.J., 1997. Modelling of runoff generation and soil moisture dynamics for hillslopes and micro-catchments. J. Hydrol. 198, 177–195. Brooks, R.H., Corey, A.T., 1964. Hydraulic Properties of Porous Media. Colorado State University. Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24, 755–769. Celia, M.A., Bouloutas, E.T., Zarba, R.L., 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26, 1483–1496. Chakravarthy, S., Osher, S., 1983. High Resolution Applications of the Osher Upwind Scheme for the Euler Equations. American Institute of Aeronautics and Astronautics. Clement, T.P., Wise, W.R., Molz, F.J., 1994. A physically based, two-dimensional, finite-difference algorithm for modeling variably saturated flow. J. Hydrol. 161, 71–90. Corradini, C., Morbidelli, R., Flammini, A., Govindaraju, R.S., 2011. A parameterized model for local infiltration in two-layered soils with a more permeable upper layer. J. Hydrol. 396, 221–232. Corradini, C., Morbidelli, R., Melone, F., 1998. On the interaction between infiltration and Hortonian runoff. J. Hydrol. 204, 52–67. DiCarlo, D.A., 2013. Stability of gravity-driven multiphase flow in porous media: 40 Years of advancements. Water Resour. Res. 49, 4531–4544. Gastó, J.M., Grifoll, J., Cohen, Y., 2002. Estimation of internodal permeabilities for numerical simulation of unsaturated flows. Water Resour. Res. 38, 62-1–62-10. Harten, A., Lax, P., van Leer, B., 1983. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev. 25, 35–61.

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032

12

H. An, S.J. Noh / Journal of Hydrology xxx (2014) xxx–xxx

Haverkamp, R., Vauclin, M., 1979. A note on estimating finite difference interblock hydraulic conductivity values for transient unsaturated flow problems. Water Resour. Res. 15, 181–187. Hu, W., Tallon, L.K., Si, B.C., 2012. Evaluation of time stability indices for soil water storage upscaling. J. Hydrol. 475, 229–241. Jones, J.E., Woodward, C.S., 2001. Newton–Krylov-multigrid solvers for large-scale, highly heterogeneous, variably saturated flow problems. Adv. Water Resour. 7, 763–774. Kelley, C.T., 1995. Iterative methods for linear and nonlinear equations. SIAM. Kurázˇ, M., Mayer, P., Lepš, M., Trpkošová, D., 2010. An adaptive time discretization of the classical and the dual porosity model of Richards’ equation. J. Comput. Appl. Math. 233, 3167–3177. Lott, P.A., Walker, H.F., Woodward, C.S., Yang, U.M., 2012. An accelerated Picard method for nonlinear systems related to variably saturated flow. Adv. Water Resour. 38, 92–101. Louaked, M., Hanich, L., 1998. TVD scheme for the shallow water equations. J. Hydraul. Res. 36, 363–378. Mingham, C.G., Causon, D.M., 1998. High-resolution finite-volume method for shallow water flows. J. Hydraul. Eng. 124, 605–614. Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12, 513–522. Paniconi, C., Putti, M., 1994. A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems. Water Resour. Res. 30, 3357–3374. Rawls, W.J., Gish, T.J., Brakensiek, D.L., 1991. Estimating soil water retention from soil physical properties and characteristics. In: Stewart, B.A. (Ed.), Advances in Soil Science. Springer, New York, pp. 213–234.

Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. J. Appl. Phys. 1, 318–333. Ridolfi, L., D’Odorico, P., Porporato, A., Rodriguez-Iturbe, I., 2003. Stochastic soil moisture dynamics along a hillslope. J. Hydrol. 272, 264–275. Roe, P.L., 1986. Characteristic-based schemes for the Euler equations. Annu. Rev. Fluid Mech. 18, 337–365. Sweby, P.K., 1984. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM J. Numer. Anal., 995–1011. Szymkiewicz, A., 2009. Approximation of internodal conductivities in numerical simulation of one-dimensional infiltration, drainage, and capillary rise in unsaturated soils. Water Resour. Res. 45, W10403. Szymkiewicz, A., Helmig, R., 2011. Comparison of conductivity averaging methods for one-dimensional unsaturated flow in layered soils. Adv. Water Resour. 34, 1012–1025. Toro, E.F., 2000. Shock-capturing Methods Shallow Flows. Wiley, Chichester. Van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892. Vereecken, H., Huisman, J.A., Bogena, H., Vanderborght, J., Vrugt, J.A., Hopmans, J.W., 2008. On the value of soil moisture measurements in vadose zone hydrology: a review. Water Resour. Res. 44, W00D06. Warrick, A.W., 1991. Numerical approximations of darcian flow through unsaturated soil. Water Resour. Res. 27, 1215–1222. Western, A.W., Grayson, R.B., Blöschl, G., 2002. Scaling of soil moisture: a hydrologic perspective. Annu. Rev. Earth Planet. Sci. 30, 149–180.

Please cite this article in press as: An, H., Noh, S.J. High-order averaging method of hydraulic conductivity for accurate soil moisture modeling. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2013.12.032