Developing a three-dimensional coupled model of pipe-matrix subsurface flow

Over the past two decades, many studies have reported the presence of soil pipes in hillslopes and their significant influence on rainfall-runoff processes. To analyze pipe flow mechanisms which have complex flow dynamics and interaction with water in the surrounding soil, this study proposed a numerical simulation model which combined a slot model with a threedimensional saturated-unsaturated subsurface flow model. Soil matrix flow and pipe flow were regarded as separate flow systems and calculated using the individual governing equations, which are Richards equation and the dynamic flow equation. To validate the model, the simulations were conducted for three different conditions (no pipe, open pipe and closed pipe) and showed good agreement with experimental observation data.


INTRODUCTION
Many researchers have reported the development of chains of interconnected macrospores called soil pipes nearly parallel to the ground surface in hillslopes.First, rainfall infiltrates into the soil matrix and the lateral flow to the pipe starts after the soil matrix becomes saturated.Then the pipe flow begins, and water can flow out from the pipe to the soil matrix once it is filled.An experimental study using a fiberscope demonstrated that both full and partially filled conditions occurred simultaneously within the same soil pipe (Terajima et al. 2000).These soil pipes play an important role related to hillslope hydrological processes and hillslope stability.The field observations conducted by Kitahara (1992) show that pipe flow can be expressed approximately using the Darcy-Weisbach equation.Further he found that the pipe flow drains out water from the hillslope quickly and it increases the slope stability.However, when the closed pipe condition occurs, it may create high pressure zones at the lower part of the pipe, which could be a cause of slope failure.
As described above, pipe flow has a significant effect on the hydrological processes of the hillslope, and the pipe flow mechanisms have complex flow dynamics and interactions with water in the surrounding soil.However, a few studies have proposed models considering detailed interaction between soil matrix flow and pipe flow at the hillslope scale.A model proposed by Kosugi et al. (2004) treated a soil pipe as a highly permeable soil layer and both the matrix flow and the pipe flow were calculated using a saturated-unsaturated flow model.Tsutsumi et al. (2005) developed a model describ-ing three-dimensional water flow in a hillslope with soil pipes.In their model, matrix flow and pipe flow are calculated using the Richards and Manning equations, respectively, while simultaneously considering the interaction between these two flow systems with iterative computations.An et al. (2007) proposed a coupled model combining a two-dimensional saturated-unsaturated flow model with a slot model which calculated pipe flow based on the dynamic flow equation.This coupled model regards matrix flow and pipe flow as separate flow systems as in the model developed by Tsutsumi et al. (2005).However it does not conduct iterative computations between the two systems.The twodimensional saturated-unsaturated flow model and slot model alternately calculate the flows in their respective domains, using each other's state variables as boundary conditions.The aim of this study is to extend the model developed by An et al. (2007) to a three-dimensional model and conduct numerical simulations to validate the extended model.

Basic Concept of the Model
The basic concept of the model developed in this study is the same as that of An et al. (2007).As the hydraulic characteristics of pipe flow and soil matrix flow are considerably different, we treat each flow as a separate flow system.Pipe flow is calculated by the slot model and soil matrix flow is calculated by the saturated-unsaturated flow model.The slot model is capable of handling both open-channel and surcharged flows with an identical flow equation and being widely used in the calculation of the urban sewerage network (Watanabe et al. 1989).The complex flow dynamics within the hillslope are then calculated by coupling the two models as follows (Figure 1): First, the saturatedunsaturated flow model calculates soil matrix flow using the water depth of pipe region as a pressure head boundary condition which is computed by the slot model.Then the interflow between the pipe and soil matrix is estimated using Darcy's law.Finally, the slot model calculates pipe flow using the estimated interflow as a lateral boundary condition.The following two sections give brief descriptions of the saturatedunsaturated flow model and the pipe flow model.

Saturated-unsaturated Flow Model
The saturated-unsaturated flow is governed by Richards equation.The three-dimensional Richards equation is written as: where is the soil moisture content, is the pressure head, and K is the hydraulic conductivity, w is the gradient of the slope, t is the time, x is parallel to the slope, z is perpendicular to the slope, and y is the horizontal coordinate.
Equation ( 1) is solved using the modified Picard method (Celia et al. 1990) and a finite volume method.The hydraulic conductivity of the boundary between adjacent control volume units is defined as (Shiraki 1998): where i, j, k is the hydraulic head of the control volume (i, j, k).

Pipe Flow Model
The pipe flow is calculated by the slot model which assumes that a pipe has a slot of which width is vanishingly small and frictional resistance is negligible (Figure 2).The fluid motion within the slot model is described by the basic equation of the open-channel flow which is written as: where h is the depth of water measured from the bottom of the pipe, g is the gravitational acceleration, Q is the discharge, A is the area of flow, w is gradient of the pipe which is the same as the gradient of the slope in this study, Sf is the frictional gradient and q is the lateral flow discharge.Sf is written as where n is Manning roughness coefficient, C is Chezy coefficient, and s is the length of wetted perimeter.In Eqs. ( 3) and ( 4), the flow between soil matrix and pipe is treated as lateral flow.Equations ( 3) and ( 4) are solved using the Preissmann four-point finite difference scheme.In this solution, the terms A/ h and s/ h arise.When h reaches zero, the terms become undefined ( A/ h h=0 , s/ h h=0 ).In order to avoid this difficulty, we define the cross sectional shape of the pipe as shown in Figure 3.There are two shapes, one is almost round but has a small notch on the bottom, and the other is square.For both of these, if h reaches zero, the values of A/ h and s/ h become constant ( A/ h h=0 const, s/ h h=0 const).

Simulation Condition
In this study, the simulations were compared with laboratory experiments conducted by Uchida et al. (1995).For comparison, the simulations were conducted using a two-dimensional model developed by An et al. (2007) and the three-dimensional model proposed in this study.The conditions of the experiment are as follows.
A flume, 70 cm long and 7 cm wide, is inclined at 15°and filled with Toyoura standard sand to a thickness of 10 cm.To simulate a soil pipe, an acrylic pipe of 1 cm outer diameter, 0.8 cm inner diameter, and 30 cm length is used.Four tiny perforations of 0.2 cm diameter are made in the walls of the acrylic pipe at 2 cm intervals and the pipe is wrapped with a cotton cloth.The upstream end of the pipe is filled with silicon and closed in order to prevent the sand particles from entering the -53 -   pipe.The water level at the lower end of the flume is fixed, and water is supplied to the upper tank at a constant rate of 0.5 cm 3 s 1 by using a rotary pump.
The relationship between the volumetric soil water content , the pressure head and the hydraulic conductivity function K( ) is expressed by the following equations (Kosugi 1994(Kosugi , 1996)): where s and r are the saturated and residual water contents respectively, m is the pressure head corresponding to the median pore radius, is a dimensionless parameter to characterize the width of the pore-size distribution, and F denotes the complementary normal distribution function, defined as: Kosugi et al. (2004) simulated the same experiments as in this study.The values of s, r, m and used in this study are from Kosugi et al. (2004) and shown in Table I.
The simulation was conducted for three different cases.In the first case, the simulation was performed without a pipe (no pipe condition).In the second case (open pipe condition), the pipe is buried so that its outlet is connected directly to the free water in the lower end of the flume (Figure 4-a).In the third case (closed pipe condition), the outlet of the pipe is placed at a 15 cm upslope from the lower end of the slope (Figure 4-b).In the case of the simulation of no pipe condition, the saturated-unsaturated flow model is used without being coupled with the slot model.In the simulations of the open and closed pipe conditions, the coupled model is used.
For the numerical stability of the simulation, a very small constant amount of discharge was added from the upstream end of the pipe.That discharge accounted for only a small fraction (< 1%) of the outflow rate at the downstream end of the pipe throughout the simulations.Also, we tested the two types of cross-sectional shape of pipe shown in Figure 3 and found that the calculation with square shape (Figure 3-b) was stable and robust.The following results are obtained using a square shape with 0.007 × 0.007 m.The grid size of x and z were 0.005 m and the time step was 10 sec in this study.The value of C shown in Table I was calibrated using the threedimensional model.For the initial condition of the simulations, we used the steady-state condition with the water level at the downstream end fixed and no water supplied from the upstream end, as in the experiment by Uchida et al. (1995).

Results of Numerical Simulations
Figure 5 shows the water surface profiles obtained from the experiment, the two-dimensional (2D) model, and the three-dimensional (3D) model.In the no pipe condition, the water surface profile of the simulation gradually increases in the upslope region with time and it matches well with the experiment.In the open pipe condition, the water surface profile of the 3D model gradually increases with time and is in reasonable agreement with that of the experiment.However, in the result of the 2D model, the profile of the lower part does not match with that of the experiment.In the closed pipe condition, the water surface profile of the 3D model is similar to that of the experiment.However, the overall water profile of the 2D model is much lower than that of the experiment and the 3D model.
These simulations show that the 2D model tends to predict the lower water table compared to the 3D model.Figure 6 is the ratio of the pipe flow discharge to the total discharge (pipe flow + matrix flow) at the lower end of the flume for the open pipe simulation.Approximately 30% of the water was discharged through the pipe in the three-dimensional simulation.On the other hand, almost all of the water was discharged through the pipe in the two-dimensional simulation.Furthermore, we repeated the 2D simulations with different values of C, but we were not able to obtain good agreement for the water table between the simulations and the experiment.These facts suggest that the effect of the pipe is overestimated in the two-dimensional model.Also, even in this small experimental flume, the flow is three-dimensional.Figure 7 shows the water flux vector and water surface profile at the y z plane.At the upper end of the open pipe, water flows into the pipe from the surrounding soil (Figure 7-a), while water flows out from the pipe to the soil at the lower end of the closed pipe (Figure 7-b).Accordingly, the twodimensional model is insufficient for dealing with the interaction between pipe flow and soil matrix flow and the three-dimensional model needs to be used to describe water dynamics in the heterogeneous soil layer.

Figure 2 .
Figure 2. Schematic representation of the slot model.

Figure 3 .
Figure3.Cross-sectional shape of the pipe in the slot model.
Figure 4. Schematic diagram of the experimental setup.(a) open pipe and (b) closed pipe.

Figure 5 .Figure 6 .
Figure 5.Comparison of the water surface profile between simulation results and experiment observations (observation interval was 10 cm).All dimensions are in cm.(a) no pipe condition, (b) open pipe condition, and (c) closed pipe condition.

Figure 7 .
Figure 7. Water flux vector and water surface profile in the y z plane (a) at the upper end of open pipe and (b) at the lower end of closed pipe.

Table I .
Parameter values.The saturated-unsaturated flow is calculated using Richards equation.The simulations were conducted for three different conditions (no pipe, open pipe, and closed pipe) using the coupled model.The simulations were compared with experimental observations and the result of the 3D model showed reasonable agreement.The 2D model tends to overestimate the effect of the pipe, but the 3D model successfully describes the water dynamics of the pipe and soil matrix.The 3D model developed in this study is capable of appropriately dealing with the interaction between pipe flow and soil matrix flow.