MODELLING TRANSIENT SEEPAGE BETWEEN PARALLEL TRENCH DRAINS

 

M.Mavroulidou & R.I.Woods,

 

Dept. of Civ. Eng., Univ. of Surrey, Guildford, Surrey, GU2 5XH, UK

Email: M.Mavroulidou@surrey.ac.uk

 

M.J. Gunn,

 

School of Construction, South Bank University, Borough Road, London, SE 10AA, UK

Email: gunnm@sbu.ac.uk

 

 

ABSTRACT

Most of the existing solutions for trench drainage systems, assume saturated soil theories. It can be anticipated however, that trench drainage involving lowering of the water table might produce desaturation of the soil. The influence of the unsaturated zone on the rate at which the transient seepage develops can be significant when shallow groundwater table movement is involved and should be considered by the trench drainage model. Due to complications involved with unsaturated hydraulic property calculations, simplified forms of these properties have often been used in the geotechnical engineering literature, for the numerical solution of flow in variably saturated porous media. The writers assess the performance of some of these simplifying assumptions, applied to the numerical modelling of trench drainage. Comparisons are made with sandbox model results found in the literature, as well as results derived by the writers' model, that incorporates possible desaturation of the materials based on data from real soils, rather than simplified forms of the unsaturated hydraulic properties. It is shown that the rate at which the transient seepage develops might be considerably underestimated when simplifying assumptions about the storage coefficient are made, and especially in coarse soils. Thus, the importance of using realistic values for the storage term variation, when modelling transient drainage conditions, is demonstrated.

 

Keywords: numerical modelling, trench drains, variably saturated soil, transient response, storage coefficient

1. INTRODUCTION

Trench drains have long been used in civil engineering to perform subsurface drainage. Their function as subsurface drainage systems is to enhance soil shear strength, which might involve reduction of pore pressures (with or without groundwater table lowering), seepage force control and protection of unstable (expansive or collapsible) subgrade soils from moisture variations due to climatic effect or shallow water table fluctuations. In particular, when transient conditions involving movement of the water table are considered, one of the main concerns of the geotechnical engineer designing a drainage system, is to obtain estimates of how the pore pressures and the water table will vary and their impact on the drain effectiveness with time. Transient phenomena such as consolidation of cohesive soils as a consequence of water table lowering due to drainage must also be taken into account in drainage system design, in order to prevent undesirable subsidence phenomena in inhabited areas. Successful modelling of drainage system design and drain efficiency with time is therefore an important topic in geotechnical engineering.

 

Due to the complexity of the problem and the difficulty in exact measurements of the factors involved (e.g. the boundary conditions and the soil properties), trench drain design has been based on semi-empirical formulae up to recently. Analytical and in the last two decades, numerical solutions have also been used to estimate flow rates and pressure distributions associated with sub-surface drainage systems. For instance, soil consolidation theories (e.g. Biot [1941]) can be used to estimate transient flow between drains, where the flow is confined (e.g. Di Maio and Viggiani [1987]). The equations obtained are similar to those found in other application areas, for example the conduction of heat through solids. Use can therefore be made of standard solutions for heat conduction (e.g. Carslaw and Jaeger, [1948]) as for example in Bromhead [1984], who applies an analytical solution for heat flow provided by Carslaw and Jaeger [1948], to the solution of a fully penetrating trench drain problem, relevant to slope stabilisation. These solutions are all linear and so the principle of superposition can be employed. Thus, to represent a potential distribution due to complex boundary conditions, the solution can be found by superposing the basic solutions, i.e. those resulting when considering the partial effect of one boundary condition at a time (e.g Bromhead [1984]). In cases where a simple reduction of the pore pressure is not sufficient, but where lowering of the water table is also necessary (e.g. for road subsoils where a shallow water table is present), the trench drainage model also has to address the issue of how the moving water table is handled. Analytical solutions to this problem are often based on simplifying assumptions such as that of predominantly horizontal flow (e.g. Maugeri and Motta [1987] solving one-dimensional transient trench drainage seepage as a sequence of steady states, using Dupuit-Forchheimer assumptions), homogeneous and isotropic soil.

 

The results of the previously mentioned solutions are often used to provide charts (e.g. giving the piezometric head reduction at any level, drain efficiency etc.), to be used for drain system design (e.g. Bromhead [1984]). Nevertheless, there are usually significant discrepancies between observed field data and predicted curves (Hutchinson [1977]). Apart from other well known difficulties (i.e. the selection of appropriate boundary conditions and the measurement of hydraulic properties of the soil), these discrepancies can also be due to the fact that these solutions inherently adopt assumptions of saturated flow. Indeed, when the soil desaturates (as can happen in trench drainage cases, involving lowering of the water table), we can expect different transient response of the system, than that based on the elastic strains of the soil skeleton assumed in consolidation theory. The influence of the unsaturated zone is anticipated to be important in the case of shallow groundwater table movement (lowering or seasonal fluctuation), often involved in trench drainage problems.

 

Unsaturated zone modelling is a relatively new topic for the geotechnical engineering domain, as opposed to other fields such as soil science and hydrology. Even when an unsaturated flow model is introduced into geotechnical numerical codes, often simplifications are made on the grounds of uncertainties considering hydraulic property variation in the unsaturated zone. For instance, Lam et al [1987] and Fredlund and Rahardjo [1993] proposed the use of a constant storage coefficient term in the unsaturated zone, in the order of 10^-3 m^-1 and 10^-2 m^-1 respectively. Moreover, linear (although sometimes logarithmically linear) hydraulic conductivity curves are used by these authors for the unsaturated zone, allowing the hydraulic conductivity to be reduced up to a certain limit (usually between 10^-4 and 10^-7 times the saturated hydraulic conductivity value). These simplified curve forms are supposed to be sufficiently satisfactory for engineering purposes. The writers have instead, a formulation considering flow in both saturated and unsaturated zones, in which the hydraulic conductivity and storage coefficient are taken as continuous functions of the pressure head, based on real soils data. The scope of the paper is to perform comparisons between the writers approach, experimental results found in the literature and a simplifying approach to the transient solution of seepage in variably saturated soils, allowing the use of a constant storage coefficient in the unsaturated zone. The purpose of this exercise is to assess the validity of the latter approach.

 

2. PROPOSED NUMERICAL MODEL

The governing equation for transient flow through an unsaturated-saturated soil system is:

(1)

where h and y are the total and pressure head respectively, K_x, K_y and K_z are the hydraulic conductivities in the x, y and z directions respectively, given as functions of the pressure head y and S is the storage coefficient. When y>0 (saturated zone), changes in soil volume are related to the elastic compressibility of the soil skeleton. Thus, when y>0 the storage coefficient S is equal to the elastic storage coefficient S_s. When y<0 the soil may become desaturated and the storage coefficient is obtained from the slope of an experimentally determined soil water characteristic curve (relating volumetric water content to pressure head). The non-linear relationship of volumetric water content and hydraulic conductivity to pressure head in the unsaturated zone, is expressed respectively by:

(2)

and

(3)

The storage coefficient for y<0 is obtained by differentiation of eqn(2) with respect to the pressure head. The parameters a₁,n₁ and a₂,n₂ involved in these relationships are back-fitted from experimental data. In the present analyses, the experimental water characteristic and hydraulic conductivity curves for a silty sand given by Vauclin et al.[1976] have been used (Figs 2(a) and 2(b)).

 

For the derivation of eqn(1) validity of Darcy's law in both the saturated and the unsaturated zones has been assumed throughout the analysis. The dependence of both the hydraulic conductivity as well as the storage coefficient on the head variations causes eqn (1) to be strongly non-linear. For the linearisation of the equation a Picard iterative scheme was adopted and a backward Euler finite-difference scheme was used for time integration.

 

3. TRENCH DRAINAGE ANALYSES

The problem solved herein, consists of a 2m high and 6m long flow domain drained by two parallel drains. The total head is initially 1.45m at both the upstream and downstream faces, when it is suddenly lowered at 0.75m at both drains. Absence of rainfall was assumed. In reality, however, boundary conditions will change continuously between the two extremes of dry surface conditions (zero flux) and ponding, depending on the succession of dry and rainy periods. The geometry and boundary conditions do not represent any actual trench drain section, but they were selected as such, to correspond to those assumed by Vauclin et al [1976], such that comparison with these authors' sandbox model results can be made. The geometry and boundary conditions of the problem are shown in Fig 1. A silty sand soil was assumed, with the properties presented in Figs.2(a)-2(c). Because of symmetry only half of the problem domain (i.e. 3m) was considered in the finite element discretisation.

 

 

Fig. 1 Geometry of the problem analysed here

 

 

Fig. 2a) Hydraulic conductivity function used here

 

Fig. 2b) Volumetric water content function used here

 

 

Fig. 2c) Storage coefficient function used here

 

Four sets of solutions are presented herein, one corresponding to the writers' model and three others presenting solutions derived by the writers' program in which constant storage coefficient assumptions were made. Namely, constant storage coefficients equal to 10 times the value of the elastic storage coefficient S_s=0.0001m^-1(similar to the storage coefficient used in Lam et al. [1987]), 100 times the value of the elastic storage coefficient (similar to storage coefficient values assumed in Fredlund & Rahardjo [1993]) and 1000 times the value of the elastic storage coefficient (storage coefficient value corresponding to the order of magnitude of the specific yield of the soil, often used to quantify the flow volumes in the unsaturated zone) were assumed. The other simplifying assumption often made, that of linear hydraulic conductivity-suction curves has been discarded herein, and the hydraulic conductivity curve presented in Fig. 2(a) (based on experimental soil data), has been used for all solutions. The purpose was therefore to focus on the impact of assuming a constant storage coefficient in the unsaturated zone.

 

The results showing the comparative phreatic surface positions at various time levels, are represented in Figs 3(a) and 3(b). From the figures it can be seen that the constant storage coefficient model results differ significantly from the writers' varying storage coefficient model. The former tend in general to speed up the rate at which transient seepage develops with respect to the results derived by the writers' varying storage coefficient model. Thus, steady state configuration (corresponding to a flat water table, of elevation equal to that of the water elevation in the drains, since infiltration rate on the soil surface is assumed to be zero), is shown to be obtained as soon as in t=0.1h (Fig. 3(a)) for a constant storage coefficient in the order of 10^-3 m^-1 and at t=0.5h for a constant storage coefficient in the order of 10^-2 m^-1, when the writers' varying coefficient model has only approached (but not yet reached) steady state at 10h. Better agreement between the two models is obtained for constant storage coefficients in the order of 10^-1m^-1 but the rate at which the transient seepage develops according to the constant coefficient model is still faster. Notice (Fig. 3(a) and 3(b)) that the writers' varying storage coefficient model results for a sandy soil, compare favourably with the sandbox model results presented in Vauclin et al. [1976].

 

 

 

Fig. 3a) and Fig. 3b) Comparative results for a sandy soil

 

The general discrepancy between the two sets of results might be explained as follows. As can be seen in Fig. 2(c), the storage coefficient in the unsaturated zone is initially increasing by several orders of magnitude with respect to the elastic storage coefficient value (equal to 10^-4 m^-1), it then presents an inflection point and subsequently starts decreasing by several orders of magnitude. A high storage coefficient term implies that the water retention capacity of the unsaturated soil zone is high. By keeping a constant, rather low storage coefficient, one fails to represent this variation of soil storage capacity with suction, accounting for high volumes of water retained in the soil. It implies, therefore, a faster release of water for a wide range of suctions, and this can be the explanation why the constant coefficient method generally approaches a steady state more rapidly than the varying coefficient method. Indeed, for the sandy soil assumed herein, a constant storage coefficient in the order of 10^-2 m^-1 might be a reasonable approximation for pressure head values up to -0.05m (see Fig. 2(c)). Subsequently, the storage coefficient gradually attains values in the order of 10^-1 m^-1, to decrease again abruptly after a pressure head value of about -0.80m. Thus, for the majority of pressure values developed between 0.1h and 10h according to the variable storage coefficient model (not shown here, due to space limitations), the sand has a high storage coefficient in the order of 10^-1 m^-1and up to about 0.7m^-1. That is the reason why, the constant coefficient in the order of 10^-1 m^-1 results do not differ from the varying coefficient results as much as the other constant coefficient results. Nevertheless, the transient movement of the water table is still faster than that of the varying coefficient model (Fig. 3(b)), probably because the storage coefficients used in the latter model are still generally higher than 0.1m^-1, i.e. the constant storage coefficient value assumed here (according to Fig. 2(c) they might be up to about seven times this value).

 

4. CONCLUSION

Unsaturated zone plays an important role in the transient response of a soil. For the successful simulation of transient drain seepage not only the hydraulic conductivity, but also the storage coefficient value will be of great importance. Simplifying assumptions often made about the value of the latter, fail to represent the complex processes developed during unsaturated flow and may lead to discrepancies with observed soil behaviour, due to inadequate representation of the drainage process. Such assumptions may speed up the rate at which the transient seepage develops, especially in coarse grained soils, leading for instance to underestimation of the time needed for the drainage system to be effective or to non-adequate representation of the seasonal variations of the shallow water table. Proper simulation of the physics of unsaturated flow gives the possibility to design drainage systems with more confidence.

 

REFERENCES

Biot, M.A.,"General theory of three dimensional consolidation", J. of Appl. Physics ,12,1941,pp.155-164.

Bromhead E.N., An analytical solution to the problem of seepage into counterfort drains, Canadian Geotechnical Journal, 21, 1984, pp. 657-662.

Carslaw, H.S. and Jaeger, J.C. Conduction of Heat in Solids, (2ed), Clarendon Press, Oxford, 1959.

Di Maio C., and Viggiani C., Influence of intermittent rainfall on effectiveness of trench drains, in Proc. of the 9^th Int. Conf. on Soil Mech. and Foundation Eng., Balkema, Roterdam, 1987, pp. 149-152

Fredlund D.G. and Rahardjo H.,Soil Mechanics for Unsaturated Soils, J. Wiley & Sons, N. York,1993.

Hutchinson J.N., Assessment of the effectiveness of correction measures in relation to geological conditions and types of slope movement, Bull. of the Int. Ass. of Eng. Geology,16,1977, pp131-155.

Lam L., Fredlund D.G., Barbour S.L., Transient seepage model for saturated-unsaturated soil systems: a geotechnical engineering approach, Canadian Geotechnical Journal, 24, 1987, pp. 565-580.

Maugeri M. and Motta E., One-dimensional modelling to assess trench drainage with time, in Proc. of the 9^th Int. Conf. on Soil Mech. and Foundation Eng., Balkema, Roterdam, 1987, pp. 189-192

Vauclin M., Khanji D., Vachaud G., Etude expérimentale et numérique du drainage et de la recharge des nappes à surface libre, avec prise en compte de la zone non saturée, Journal de Mécanique,15(2), 1976, pp.307-348.