Repositorio Institucional de Documentos

Numerical simulation of one-dimensional transient vertical flow in variably saturated soils Daniel Caviedes Voulli`eme M´aster en Mec´anica Aplicada Programa Oficial de Posgrado en Ingenier´ıa Mec´anica y de Materiales September 2010 Supervisor: Dr. Pilar Garc´ıa Navarro Term 2009-2010 Centro Polit´ ecnico Superior Universidad de Zaragoza Abstract iii Numerical simulation of one-dimensional transient vertical flow in variably saturated soils Abstract Water flow in variably saturated (saturated/unsaturated) soils is commonly modeled by means of Richards’ equation. This equation has no general analytical solution and the use of numerical approximations is necessary. It can be presented in three physically equivalent forms which are based on different variables and show different mathematical properties. In this work, these forms are derived from the general mathematical model and analyzed from a numerical perspective in order to understand the interactions between the differential equations and the numerical methods required in each case. The goals of this work are, on the one hand, to describe the physical and mathematical reasoning which leads to the formulation of the general mathematical model of flow in porous media, the discussion of the concepts and assumptions which allow to develop Richards’ equation, and on the other hand, to establish the properties and limitations of several numerical schemes to approximate the solutions of flows in variably saturated soils. The approach for the mathematical model is to average a microscopic, singlephase flow equation into a macroscopic scale which allows to describe porous media in a practical way, and to consider the necessary assumptions to state Richards’ equation as a particular flow case. The porous media constitutive model completes the mathematical model. In this work the Mualem-van Genuchten model and some variants are included. For the numerical model, several schemes are developed for the 1D Richards’ equation in the vertical direction. Explicit and implicit centered finite difference schemes are used in this work. The key numerical aspects of interest are those of mass conservation, stability and efficiency. Another key aspect, which is not only numerical is that of continuity from unsaturated into saturated regimes. The constitutive models affect the numerical schemes and some issues arise because of the high non-linearity of the functions, in particular the hydraulic conductivity function. Appropriate discretization of hydraulic conductivity for estimation of flux between numerical cells is a sensible issue wich has been studied by many authors and is treated in this work. All of these issues are analyzed individually and as interrelated problems in the schemes. Validation and test cases are presented and the response of the model to different problems and parameters is examined. From them, it is concluded that the explicit and the implicit schemes based on the mixed form of Richards’ equation are better suited for unsaturated problems. For variably saturated problems, the implicit scheme based on the mixed form is the best choice, since the explicit model cannot solve saturation conditions. Conditional stability of the explicit model affects negatively its performance in certain cases, which also leads to the conclusion that the implicit scheme is more efficient and realiable. Resumen v Simulaci´on num´erica unidimensional de flujos transitorios verticales en suelos con saturaci´on variable Resumen El flujo de agua en suelos con saturaci´on variable (saturado/no saturado) es comunmente modelizado por medio de la ecuaci´on de Richards. Dicha ecuaci´on no tiene una soluci´on an´alitica general, y por tanto es necesario el uso de aproximaciones num´ericas. La ecuaci´on puede presentarse en tres formas, las cuales son f´ısicamente equivalentes, pero basadas en distintas variables, y que muestra comportamientos matem´aticos distintos. En este trabajo, dichas formas se obtienen a partir del modelo matem´atico general, y se analizan desde la perspectiva num´erica con el fin de comprender las interacciones entre las formas de la ecuaci´on diferencial y los m´etodos num´ericos aplicables en cada caso. Los objetivos de este trabajo son, por una parte, describir el razonamiento f´ısico y matem´atico que lleva a la formulaci´on del modelo matem´atico general de flujo en medios porosos, la discusi´on de los conceptos y supuestos que permiten formular la ecuaci´on de Richards, y por otra, estudiar las propiedades y la aplicabilidad de los m´etodos num´ericos para su soluci´on. El enfoque utilizado para el modelo matem´atico es el de promediar una ecuaci´on microsc´opica de una sola fase, a una escala macrosc´opica que permita describir el medio poroso de una forma pr´actica y, posteriormente, considerar los supuestos que permiten formular la ecuaci´on de Richards como un caso particular de flujo. El modelo constitutivo del medio poroso completa dicho modelo matem´atico. En este trabajo se utiliza el modelo de Mualem-van Genuchten as´ı como una de sus variantes. Para el modelo num´erico, varios esquemas num´ericos fueron formulados para la ecuaci´on de Richards unidimensional, en la direcci´on vertical. En este trabajo se utilizan esquemas expl´ıcitos e implicitos con diferencias finitas centradas. Los aspectos de inter´es desde la perspectiva num´erica son la conservaci´on de masa, estabilidad y eficiencia computacional. Adicionalmente, un tema que no es ´unicamente num´erico es la continuidad y aplicabilidad tanto en la regi´on saturada como en la regi´on parcialmente saturada. Los modelos constitutivos de suelo llevan al modelo num´erico una serie de dificultades por las caracter´ısticas de las funciones no lineales, en particular la conductividad hidr´aulica. La discretizaci´on cuidadosa the la conductividad entre las celdas es un tema de importancia, el cual ha sido estudiado por muchos autores y se incluye tambi´en en este trabajo. Todos estos aspectos se analizan en s´ı mismos y como problemas interrelacionados dentro de los esquemas num´ericos. Se presentan pruebas de validaci´on y casos test y se examina la respuesta de los modelos a distintos problemas y par´ametros. De dichas pruebas se puede concluir que los esquemas expl´ıcitos e impl´ıcitos basados en la forma mixta de la ecuaci´on de Richards son m´as apropiados para la soluci´on de problemas no saturados. Para condiciones de saturaci´on variable, el esquema impl´ıcito basado en la forma mixta es la mejor opci´on, dado que el modelo expl´ıcito no es capaz de resolver condiciones de saturaci´on. La estabilidad condicionada del modelo expl´ıcito tambi´en afecta de forma negativa a su eficiencia computacional en algunos casos, lo cual apoya la conclusi´on general de que el esquema impl´ıcito es m´as robusto y eficiente. Contents Abstract iii Resumen v Introduction 2 1 Mathematical model and governing equations 4 1.1 Microscopic equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Averaging rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.1 Averaging definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Average of the time derivative . . . . . . . . . . . . . . . . . . . . . . 6 1.2.3 Average of the spatial derivative . . . . . . . . . . . . . . . . . . . . . 7 1.3 Macroscopic equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Macroscopic mass balance equations . . . . . . . . . . . . . . . . . . . . . . . 8 1.5 Mass balance in a non deformable, variably saturated porous medium . . . . 9 1.6 Macroscopic momentum equation . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.7 Richards’ Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.8 Unsaturated soil constitutive model . . . . . . . . . . . . . . . . . . . . . . . . 14 1.8.1 Mualem-van Genuchten Model . . . . . . . . . . . . . . . . . . . . . . 15 1.8.2 Modified Mualem-van Genuchten Model . . . . . . . . . . . . . . . . . 15 2 Numerical model 18 2.1 Spatial and temporal discretization . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2 Explicit formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2.1 Mixed scheme (EMC) . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.2 Pressure-based scheme (EP) . . . . . . . . . . . . . . . . . . . . . . . . 19 viii Contents 2.2.3 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3 Implicit formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.1 IP Scheme .................................. 20 2.3.2 IMC scheme ................................. 21 2.3.3 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3.4 Convergence and under-relaxation . . . . . . . . . . . . . . . . . . . . 24 2.4 Scheme properties ................................. 25 2.4.1 Solution method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.4.2 Transition from unsaturated to saturated . . . . . . . . . . . . . . . . 25 2.4.3 Mass conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.4.4 Stability ................................... 27 2.4.5 Efficiency .................................. 28 2.5 Computation of intercell conductivity Ki±1/2. . . . . . . . . . . . . . . . . . 29 2.6 Mass Balance Error Assesment . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3 Validation and test cases 32 3.1 Warrick’s Analytical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2 Test cases ...................................... 37 3.2.1 Test Case 1: Impervious boundaries . . . . . . . . . . . . . . . . . . . 37 3.2.2 Test Case 2: Downward saturation in semi-infinite soil . . . . . . . . . 38 3.2.3 Test Case 3: Downward saturation with water table . . . . . . . . . . 39 3.2.4 Test Case 4: Downward partial saturation process . . . . . . . . . . . 40 3.2.5 Test Case 5: Downward full saturation process . . . . . . . . . . . . . 43 3.2.6 Test Case 6: Downward drying process with water table . . . . . . . . 43 3.2.7 Test Case 7: Downward drying process in semi-infinite soil . . . . . . 44 4 Conclusions and further research 46 4.1 Conclusions ..................................... 46 4.2 Further research .................................. 47 Bibliography 48 A Numerical schemes formulations 52 A.1 Explicit Mixed Scheme ............................... 52 Contents ix A.2 Explicit pressure based scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 52 A.3 Implicit Presure based scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 A.4 Implicit Mixed Conservative scheme . . . . . . . . . . . . . . . . . . . . . . . 54 B Stability Analysis 56 B.1 EMC Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 B.2 IMC Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 1.2. Averaging rules 5 Thus, the rate of change of ein the neighborhood of a point is given by ∂eα ∂t =−∇(eαVα+jα) + ρΓE α(1.1) 1.2 Averaging rules In order to transform equation (1.1) into a macroscopic equation which is valid in a volume Uowhich contains several phases, it is necessary to introduce some averaging definitions, to transform mathematical point properties into representative properties in a control volume. Consider two phases and volume Uo=Uoα +Uoβ as in figure 1.1 for the following definitions. 1.2.1 Averaging definitions Let eαbe the volumetric phase average of eαin volume Uo. eα=1 UoZ Uoα eαdUα(1.2) Let eααbe the volumetric intrinsic phase average of density eαin the α-phase fraction Uoα of volume Uo. eαα=1 Uoα Z Uoα eαdUα(1.3) By defining the volume fraction θα=Uoα Uo , both averages can be related: eα=θαeαα(1.4) Let ˚eαbe the deviation of eαof a mathematical point from the intrinsic phase average: ˚eα=eα−eαα(1.5) Let z} eαβ be the average of eover a Sαβ surface: z} eαβ =1 Sαβ Z Sαβ e dS (1.6) Let Σαβ be the specific area of Sαβ, this is Σαβ =Sαβ Uo, so that z} eαβΣαβ =1 UoZ Sαβ e dS (1.7) A relevant property of the intrinsic phase average is that it is a linear operator. Then, any quantity Gsatisfies, G1+G2 α=G1 α+G2 α(1.8) The intrinsic phase average of a product of Gdefined, following (1.3) by G1G2 α=1 Uoα Z Uoα G1G2dUα 6Chapter 1. Mathematical model and governing equations Because of (1.5), and considering that the intrinsic average of the fluctuations is zero, G1G2 α=G1 αG2 α+˚ G1˚ G2 α(1.9) 1.2.2 Average of the time derivative Consider Reynolds theorem for the extensive property E, over a volume Uoα contained by a surface Sα=Sαα +Sαβ containing such volume with normal outwards pointing vector ˆn.Soα can be considered as the sum of the contact surface between phases αand βand the surface which separates phase αwithin the volume Uoand the outside of such volume (Sαα). DE Dt Z U(t) e dU =Z U(t) ∂e ∂t dU +Z S(t) eVEˆn dS where DEG Dt =∂G ∂t +VE∇Gis the material derivative with respect to an observer moving with S(t). Assumption 1.2.1. Uoα is assumed a material volume with respect to E, hence Soα is a material surface which implies that VE=ufor surface Sαβ where uis the velocity at which Sαβ is being displaced. DE Dt Z Uoα(t) e dU =Z Uoα(t) ∂e ∂t dU +Z Sαβ (t) euˆn dS +Z Sαα(t) eVEˆn dS (1.10) On the other hand, by considering the entire volume Uo, it is possible to express the material rate of change of the extensive quantity Ewhich only exists within phase α. To do this, consider the characteristic function γαfor phase α: γα=(1 , for points within Uoα 0 , for points outside Uoα (1.11) Then, the material rate of change of Ewithin the entire volume is DE Dt Z Uo γαe dU =∂ ∂t Z Uo γαe dU +Z So γαeVEˆn dS which, when evaluating γαyields DE Dt Z Uoα e dU =∂ ∂t Z Uoα e dU +Z Sαα eVEˆn dS (1.12) By equating (1.10) and (1.12) considering instantaneously all time dependent terms, ∂ ∂t Z Uoα e dU |{z } 1  =Z Uoα ∂e ∂t dU |{z } 2  +Z Sαβ euˆn dS Note that terms 1 and 2 may be substituted by using equation (1.3) ∂ ∂t (eααUoα) = Uoα ∂eα ∂t α +Z Sαβ euˆn dS 1.3. Macroscopic equation 7 By dividing the entire equation by Uo(which is constant in time, hence can go into the derivatives), and recalling that θα=Uoα Uo and rearranging, yields θα ∂eα ∂t α =∂ ∂t (eααθα)−z} euˆn αβ Σαβ (1.13) This equation relates the average of a time derivative of eαto the time derivative of the average of eα. 1.2.3 Average of the spatial derivative From Gauss’ theorem, ZUoα ∇G dU =ZSαα Gˆn dS +ZSαβ Gˆn dS (1.14) Consider the integral over Sαα, by means of the characteristic function γαshown in (1.11) Z Sαα Gˆn dS =Z So Gγαˆn dS Using Gauss’ theroem once again Z Sαα Gˆn dS =Z Uo ∇(Gγα)dU Because Uodoes not change in space, the order of integration and differentiation can be exchanged, and furthermore, evaluating γα, Z Sαα Gˆn dS =∇Z Uo GγαdU =∇Z Uoα G dU Substituting in (1.14), Z Uoα ∇G dU =∇Z Uoα G dU +Z Sαβ Gˆn dS By means of (1.3) and (1.7) Uoα∇Gα=∇UoαGα+Uoz} Gˆn αβ Σαβ Dividing the entire equation by Uo, and because θα=Uoα Uo finally yields θ∇Gα=∇θGα+z} Gˆn αβ Σαβ (1.15) 1.3 Macroscopic equation Integrating equation (1.1) over the volume of phase α, i.e. Uoα, and dividing by the porous medium representative volume Uo, yields 1 UoZ Uoα ∂e ∂t =−1 UoZ Uoα ∇(eV+j) + 1 UoZ Uoα ρΓE 8Chapter 1. Mathematical model and governing equations By using equations (1.13) and (1.15), ∂(θeα) ∂t =−∇ θeαVα |{z} Advective flux +˚e˚ Vα |{z} Dispersive flux +jα |{z} Diffusive flux −z} [e(V−u) + j)]ˆn αβ Σαβ |{z } Surface flux +θρΓEα |{z } Source/Sink (1.16) which is the general macroscopic conservation equation for any property e. Note that the equation includes advective, dispersive and diffusive fluxes, as well as flux through the interphase surface and a source/sink term. 1.4 Macroscopic mass balance equations Consider mass mas the extensive property Ein α-phase of a single component. Hence, the intensive property ebecomes mass denstity ρ. By taking equation (1.16), and ΓE= Γm= 0 since mass is not generated within the volume, the general mass equation is obtained. ∂(θρα) ∂t =−∇hθραVα+ ˚ρ˚ Vα+jαi−z} [ρ(V−u) + j]ˆn αβ Σαβ (1.17) By applying the intrinsic phase average of a product defined in equation (1.9) to the definition of the tensor quantity jα, together with the linear operator property shown in (1.8), it is possible to obtain after some manipulation, ∂(θρα) ∂t =−∇hθραVα+ ˚ρ˚ Vαi−z} ρ(V−u)ˆn αβ Σαβ (1.18) Assumption 1.4.1. There is no mass exchange between phases αand β. Hence, surface Sαβ is a material surface respect to mass, which implies V−u= 0. ∂(θρα) ∂t =−∇hθραVα+ ˚ρ˚ Vαi (1.19) Assumption 1.4.2. For each fluid phase, the sum of the dispersive and diffusive fluxes of the total mass is much smaller than the advective term. Considering assumptions 1.4.1 and 1.4.2, equation (1.19) reduces to ∂(θρα) ∂t =−∇θραVα=−∇(ραqα) (1.20) where qα=θVαis the specific discharge of α-phase. Equation (1.20) is a mass conservation equation with predominant advection and immiscible phases. 1.5. Mass balance in a non deformable, variably saturated porous medium 9 1.5 Mass balance in a non deformable, variably saturated porous medium Figure 1.2: Three-phase representative volume Consider a REV such as the one shown in figure 1.2. Let Uocontain three phases only, one of which is a solid phase (s). Consider the other two as fluid phases: a wetting phase (w) and a non-wetting phase (n). A simple case of this is to imagine liquid water and air in soils. Let Sαbe saturation Sα=θα ηwhere ηis porosity, and consider that the specific discharge is qα=θαVα= qrα +θαVs, where qrα is the specific discharge of αphase relative to the (in the general case) moving solid with velocity Vs. Hence, the mass balance equation (1.20) of the wetting phase α=wcan be written as Swρw ∂η ∂t +ηρw ∂Sw ∂t +ηSw ∂ρw ∂t =−Swqrw∇ρw−ρwqrw∇Sw−Swρw∇qrw (1.21) In a similar way, the mass balance equation for the nonwetting phase is obtained by making α=n. Snρn ∂η ∂t +ηρn ∂Sn ∂t +ηSn ∂ρn ∂t =−Snqrn∇ρn−ρnqrn∇Sn−Snρn∇qrn (1.22) The equation for the solid phase is 1 1−η Ds Dt(1 −η) + 1 ρs Dsρs Dt =−∇Vs(1.23) where the material derivative DαG Dt =∂G ∂t +Vα∇Grefers to velocity Vα. By combining equation (1.21) and (1.22) with (1.23), the wetting phase mass balance equation is obtained, ηSw ρw Dwρw Dt +Sw 1−η Dsη Dt +ηDsSw Dt −ηSw ρs Dsρs Dt =−∇qrw (1.24) and the nonwetting phase ηSw ρn Dnρn Dt +Sn 1−η Dsη Dt +ηDsSn Dt −ηSw ρs Dsρs Dt =−∇qrn (1.25) Assumption 1.5.1. At the microscopic level the solid phase microscopic volume dU remains constant, hence the solid phase can be considered incompressible. Ddms Dt := D(ρsdUs) Dt = 0 z}| { ρs D(dUs) Dt +dUs Dρs Dt = 0 From this, it is clear that Dρs Dt = 0. 10 Chapter 1. Mathematical model and governing equations By considering assumption 1.5.1, equations (1.24) and (1.25) can be reduced to ηSw ρw Dwρw Dt +Sw 1−η Dsη Dt +ηDsSw Dt =−∇qrw (1.26) for the wetting phase, and ηSw ρn Dnρn Dt +Sn 1−η Dsη Dt +ηDsSn Dt =−∇qrn (1.27) for the non-wetting phase. For practical reasons, it is convenient to write all material derivatives relative to the solid phase. Multiplying equation (1.26) and (1.27) by ρw, together with the definition of relative specific discharge ηSw Dsρw Dt +Swρw 1−η Dsη Dt +ηρw DsSw Dt =−∇(ρwqrw) (1.28) ηSn Dsρn Dt +Snρn 1−η Dsη Dt +ηρn DsSn Dt =−∇(ρnqrn) (1.29) Assumption 1.5.2. The non-wetting phase has a constant and uniform pressure such that pn= 0. Because equations (1.28) and (1.29) are a coupled system which is related by pressures pnand pw. The system is decoupled and one of the equations can be neglected by this assumption. Assumption 1.5.3. At the macroscopic level, the solid matrix is immobile, hence Vs= 0. ∂η ∂t = 0. Hence, equation (1.23) is reduced to DsG Dt =∂G ∂t +Vs∇G |{z } =0 ⇒DsG Dt =∂G ∂t By assumption 1.5.2 the mass conservation equation for the non-wetting phase has been neglected and only the wetting-phase remains of interest. Finally, equation (1.28) is reduced to ηSw ∂ρw ∂t +ηρw ∂Sw ∂t =−∇(ρwqrw) (1.30) which describes flow of a wetting-phase within a non-deformable porous medium in partially saturated conditions, where the wetting-phase is immiscible with the non-wetting phase which is assumed at constant pressure and no sinks/sources are considered. Saturated conditions imply Sw= 1 and ∂Sw ∂t = 0, which results in the equation for saturated flow in non-deformable porous media, with the same restrictions as for equation (1.30), η∂ρw ∂t =−∇(ρwqrw) (1.31) 1.6 Macroscopic momentum equation Consider equation (1.16) for momentum, hence E=mVand e=ρV. By decomposing the total momentum flux in terms of the momentum flux relative to velocity: j=ρV+jM. Note 1.6. Macroscopic momentum equation 11 that the momentum flux relative to mass velocity jMis actually stress, hence jM=−σ. Furthermore, momentum generation is ΓM=Fwhere Fis the external body force acting on the phase. ∂θρVα ∂t =−∇"θρVαVα+˚ (ρV)˚ Vα#−σα−z} [ρV(V−u) + σ)]ˆn αβ Σαβ +θρFα(1.32) Using the identity (which is derived from the definitions of deviation and averages) ˚ (ρV)˚ Vα=Vα˚ρ˚ Vα+ρα˚ V˚ Vα+˚ρ˚ V˚ Vα on the right side of equation (1.32) and applying equation (1.9) to the left side, yields ∂θραVα ∂t + ∂θ˚ρ˚ Vα ∂t =−∇"θρVαVα+Vα˚ρ˚ Vα+ρα˚ V˚ Vα+˚ρ˚ V˚ Vα−σα# −z} [ρV(V−u) + σ)]ˆn αβ Σαβ +θρFα By combining with the mass conservation equation (1.18), manipulating somewhat and in indicial notation, θρα∂Vi α ∂t + ∂θ˚ρ˚ Vi α ∂t =−∇hθ˚ρ˚ Vi αVj α+ρα˚ Vi˚ Vj α+˚ρ˚ Vi˚ Vj αi+∇(θσijα) + θρFi α −θ∂Vi α ∂xjραVj α+ ˚ρ˚ Vj α−z} hρ˚ Vi(Vi−uj)−σijiˆnj αβ Σαβ Assumption 1.6.1. ρis constant Assumption 1.6.2. Sαβ is a material surface with respect to α-phase Assumption 1.6.3. Flow is macroscopically uniform: ∂Vi α ∂xj = 0 Considering these assumptions, allows for θρα∂Vα ∂t =−ρα∇θ˚ V˚ Vα+∇(θσα) + θρFα+z} σˆn αβ Σαβ Assumption 1.6.4. Dispersive mass fluxes are much smaller than advective mass fluxes |˚ρ˚ Vα| ≪ |ραVα|, thus may be considered as negligible. Assumption 1.6.5. Dispersive momentum fluxes are much smaller than advective momentum fluxes |˚ (ρV)˚ Vα| ≪ |ρVαVα| ≈ |ραVαVα|. With the aforementioned assumptions it is possible to write the macroscopic momentum balance equation: θραα∂Vα α ∂t +θραVα∇Vα=∇(θσα) + z } σˆn αβ Σαβ +θ ραFα 12 Chapter 1. Mathematical model and governing equations which because of (1.15) is θραα∂Vi α ∂t +θραVj α∂Vi α ∂xj =θ∂σij ∂xj α +θραFi α(1.33) At this point, it is necessary to clarify that, although variably saturated flow in porous media is in fact a multiphase system, the following reasoning is not one of a “true” multiphase system, since water flow will be thought of as uncoupled from air flow, in accordance to assumption 1.5.2. The analysis can be done as in fully coupled multiphase system, but it is unnecesary for the intended model, hence it may be approximated by a saturated analysis. Assuming microscopical isochoric motion, evaluating the stress tensor, introducing the no-slip condition in the solid-fluid surface and rearranging [5], ρh∂qi ∂t |{z} 1  +∂ ∂xjqiqj η |{z } 2 i=−η∂p ∂xj +ρg∂x3 ∂xiT∗ ji +µ∂2qri ∂x2 j |{z } 3  −−µαij Cf ∆2 f qm rj |{z } 4  (1.34) Where pis fluid pressure, qis specific discharge, ηis porosity, µis viscosity, T∗ ij and αij are tensorial properties fo the configuration of the solid-fluid surface when saturated with the phase of interest, Cfis a macroscopic dimensionless shape factor and ∆fis the ratio of void space volume to interface surface area. Term 1 is the temporal variation, 2 describes inertial forces, term 3 expresses the viscous forces due to shear inside the fluid and term 4 expresses the drag at the solid-fluid surfaces. Equation (1.34) can be transformed into a dimensionless form [5], from which the Reynolds (Re), Darcy (Da) and Strouhal (St) numbers for porous media can be formulated. The Reynolds number is a dimensionless ratio of inertial and viscous forces. Darcy’s law is considered valid when Re ≤1−10 which is usually true in groundwater flows [6]. Darcy’s number is a ratio of characteristic permeability, and the characteristic travel paths and distance. The Darcy number allows to estimate the magnitude of viscous resistance within the fluid, which is also related with the Reynolds number. The Strouhal number is the ratio between a characteristic travel time required to encounter a significant spatial change in velocity and a characteristic time required to encounter the same change in velocity (in time) in a mathematical point. In a sense, it can be interpreted as a ratio of local and convective accelerations. Consider no momentum transfer between fluid phases and considering inertial effects and resistance of flow from viscous shear inside each fluid neglegible with respect to the shear produced with the solid. In equation (1.34) this means 2  ≪ 4 and 3  ≪ 4 . This conditions are true when Re Da1 2≪1. Furthermore when the Strouhal number is small, St ≤1, term 1 may be neglected too. Then equation 1.34 reduces to qm rj =−1 µ ∆2 fη Cf (αji)−1T∗ il |{z } kjl ∂pf ∂xj +ρfg∂x3 ∂xi Assuming that tensors (αij)−1and T∗ il have the same principal directions, and assuming cartesian coordinates, yields qrα =−kα µα (∇pα+ραg∇z) (1.35) 1.7. Richards’ Equation 13 where kαis the effective permeability to the α-phase, and is a function of saturation. For the case of unsaturated water flow in soils, the air phase is usually considered with constant atmospheric pressure in the entire domain, hence an air phase equation in the manner of equation (1.35) is unnecesary. Only an equation for the wetting phase (i.e., water) is required. Assumption 1.6.6. Water viscosity and density remains (approximately) constant. For the water flow in soil, the effective hydraulic conductivity can be defined as Kw=kwgρw µw Piezometric head of water can be defined as hw=pw gρw . By using such definitions, equation (1.35) for water specific discharge in a soil may be written as qrw =−Kw∇(hw+z) (1.36) Note that hydraulic conductivity Kderives from permeability kwhich was defined from a saturated, single fluid phase approach. In section 1.8 models for Kwill be considered which depend on such saturated approach and modified by factors associated to variable saturation. 1.7 Richards’ Equation Combining equation (1.30), which describes mass conservation in a non-deformable porous medium, and (1.36) which describes momentum conservation, leads to ηSw ∂ρw ∂t +ηρw ∂sSw ∂t =∇hρwKw∇(hw+z)i Considering an incompressible wetting phase, such as water, results in Richards’ equation [30] (phase subindices have been dropped for simplicity in notation), η∂S ∂t =∇hK∇(h+z)i Considering the definition of saturation, it is possible to write Richards’ equation in terms of the wetting phase fraction θ[L3/L3], which in common groundwater terminology is called volumetric water content. ∂θ ∂t =∇K∇h+z(1.37) This equation relates the changes in water content with the primary driving forces: gravitational potential and pressure gradient and the properties of the porous medium, by means of its conductivity. From the mathematical point of view, Richards’ equation is a parabolic equation in unsaturated regime, and an elliptic equation in saturated regime [18]. Richards’ equation may be written in three forms depending on the choice of variables. These are shown here in 1D form. They are the water content form in terms of water content θ, conductivity K(θ) and diffusivity D(θ) ∂θ ∂t =∂ ∂z D(θ)∂θ ∂z +∂K(θ) ∂z (1.38) 14 Chapter 1. Mathematical model and governing equations the pressure or matric form in terms of pressure h, hydraulic capacity C(h) and conductivity K(h) C(h)∂h ∂t =∂ ∂z "K(h)∂h ∂z + 1#(1.39) and the mixed form in terms of pressure h, water content θand conductivity K(h). ∂θ(h) ∂t =∂ ∂z "K(h)∂h ∂z + 1#(1.40) where C=∂θ ∂h [1/L] is the hydraulic capacity of the soil and D=K C[L2/T] the diffusivity. This definitions allow to relate all three forms. Although it is implicit in the equations, it is worth to observe that specific discharge qwcan be thought of as the darcian velocity or flux J, which is positive upwards because of the adopted sign convention. The three forms (1.38) (1.39) (1.40) have different properties. The water content form is the conservative form, in the sense that the variable of interest is a conserved variable. Because of this it shows very good conservation properties [13]. However, water content varies only when in the unsaturated region, hence the equation is useless when flow occurs in a saturated regime. From a conceptual perspective, these equations do not show explicitly the driving forces of flow, since it is the pressure gradient which generates flow, which is associated to different water contents by the soil constitutive model. The pressure form involves only changes in pressure. Although, when coupled with the soil constitutive model, it also relates to water content. Because pressure is a continuous function from negative pressure (suction, matric potential) in an unsaturated regime to positive pressures in a saturated regime, the transiton is well handled by solving pressure. On the other hand, this equation is not written in terms of a conserved variable and, when solved by numerical methods, may show problems in conservation [13]. The mixed form relates the change in water content to the pressure gradients. It is, in the conceptual sense, better to understand the driving forces of flow. Furthermore, because the conserved variable is present in the equation, conservation is better handled with this equation [13] [29]. As written in equation (1.40) it is not continuous into the saturated region, since water content becomes constant. Note that from the physical perspective, all three forms model the same phenomenon. Mathematically, from the differential point of view conservation is not an issue, but from the numerical perspective it is, as will be discussed in the following chapter. Because in this work variably saturated soils are of primary interest the water content form is not used. It is also important to observe that all forms depend strongly on K(h) or K(θ), C(h), D(h) or D(θ) functions defined by a constitutive model for the porous media, which is essential. 1.8 Unsaturated soil constitutive model Richards’ equation in any form requires the hydraulic conductivity K(h) function and the water content θ(h) function to be known. These functions interrelate pressure, water content, conductivity and other soil properties. There are several models [34] that feed upon soil parameters to generate mathematical relations for the functions K(h), θ(h) and its derivative, 2.3. Implicit formulation 21 where the coefficients are ai=−∆t δz2Kn+1,m i−1/2(2.10) bi=Cn+1,m i+∆t δz2Kn+1,m i−1/2+Kn+1,m i+1/2(2.11) ci=−∆t δz2Kn+1,m i+1/2(2.12) fi=Cn+1,m ihn i+∆t δz Kn+1,m i+1/2−Kn+1,m i−1/2(2.13) 2.3.2 IMC scheme From the discretized mixed-form (2.8), but actually solving for pressure the Implicit Mixed Conservative (IMC) scheme is formulated. Solving for pressure allows for the scheme to transition from unsaturated to saturated regimes, hence it is actually a solution for variably saturated flows. The highly non-linear nature of Cis a key factor in mass conservation when solving pressure, which is not a conserved variable. Although it is correct to define C=∂θ ∂h when dealing with continous functions, it is necessary to be carefull when evaluating them in a discrete way. The approximation ∂θ ∂t ≈θn+1,m+1 −θn ∆t≈Cn+1,m hn+1,m+1 −hn ∆t disregards the fact that Cis a function of hand hence a function of t, which in the discrete form is not an accurate approximation. The time derivative of θshould consider the chain rule for the time derivative of h, both in time stepping form nto n+ 1 and from mto m+ 1. If this derivative is not considered the scheme shows poor mass conservation, as clearly shown by Celia et al. [13]. In order to consider this, Celia et al. showed that the use of the first order Taylor polynomial applied to the time derivative of θ, around hn+1,m is a good solution, for the method becomes perfectly mass conservative [13]. θn+1,m+1 =θn+1,m +Cn+1,mhn+1,m+1 −hn+1,m(2.14) This results in a better approximation ∂θ ∂t ≈θn+1,m+1 −θn ∆t≈θn+1,m +Cn+1,m hn+1,m+1 −hn+1,m−θn ∆t where the derivatives in the miteration are properly considered. Writing the equation considering the Taylor polynomial and the Picard iteration results in the Implicit Mixed Conservative (IMC) scheme, aihn+1,m+1 i−1+bihn+1,m+1 i+cihn+1,m+1 i+1 =fi(2.15) Where coefficients a, b, c and the term fare ai=−∆t δz2Kn+1,m i−1/2(2.16) 22 Chapter 2. Numerical model bi=Cn+1,m i+∆t δz2Kn+1,m i−1/2+Kn+1,m i+1/2(2.17) ci=−∆t δz2Kn+1,m i+1/2(2.18) fi=Cn+1,m ihn+1,m i+∆t δz Kn+1,m i+1/2−Kn+1,m i−1/2+θn i−θn+1,m i(2.19) This scheme was first proposed by Celia et al. [13] both for finite difference and finite element schemes. It is well-known and frequently used. 2.3.3 Boundary conditions Dirichlet boundary conditions are treated as imposed pressure head conditions, while Neumman conditions are pressure gradients. Physically, imposing positive pressure head in the upper boundary can represent surface water height. Negative upper pressure head seems less natural. Imposed pressure head at the lower boundary is somewhat difficult to imagine, as it appears artificial to have pressure below the soil column which does not depend on the soil column. Neumman conditions allow to simulate a semi-infinite stratum in the soil, which responds only to the state of the column above it. In order to impose flow or zero-flow (impervious) conditions, Richards’ equation can be written in terms of flux, which allows to write the derivative of the flux, and evaluate fluxes in the intercell boundaries and to impose one of them while expressing the other in terms of the pressure gradient and conductivity as in the schemes. Flow conditions have a very clear physical meaning, even when the imposed flow is zero, which can represent impervious strata, or no infiltration from the upper boundary. Imposed pressure head When imposing pressure head at the upper and lower boundary conditions it is necessary to eliminate the first (in the case of i= 1, the lower boundary ) or last (i=N, the upper boundary) from the system of equations or matrix system, since hn+1 ior hn+1 Nrespectively, would not be unknowns, but would be part of term f. Then, it is only necessary to impose hn+1 i=hlow(t) or hn+1 N=hupp(t). Lower boundary •IP f2=Cn+1,m 2hn 2+∆t δz Kn+1,m 2+1/2−Kn+1,m 2−1/2−a2hn+1 1(2.20) •IMC f2=Cn+1,m 2hn+1,m 2+∆t δz Kn+1,m 2+1/2−Kn+1,m 2−1/2+θn 2−θn+1,m 2−a2hn+1 1(2.21) 2.3. Implicit formulation 23 Upper boundary •IP fN−1=Cn+1,m N−1hn N−1+∆t δz Kn+1,m (N−1)+1/2−Kn+1,m (N−1)−1/2−cN−2hn+1 N(2.22) •IMC fN−1=Cn+1,m N−1hn+1,m N−1+∆t δz Kn+1,m (N−1)+1/2−Kn+1,m (N−1)−1/2+θn N−2−θn+1,m N−2−cN−2hn+1 N (2.23) Imposed flow Flow at the boundaries is imposed at the lower face of cell i= 1 or the upper face of cell i=N. The first is denoted Joand the latter J∗, both positive upwards as shown in figure 2.1. The corresponding coefficients for the first or last row of the matrix equation are as shown below. Note that there is no a1term as there is no hn+1 0unknown at the lower boundary, and at the upper boundary there is no cNterm as there is no hn+1 N+1 unknown. Lower boundary b1=Cn+1,m 1+∆t δz2Kn+1,m 1+1/2(2.24) c1=−∆t δz2Kn+1,m 1+1/2(2.25) •IP f1=Cn+1,m 1hn 1+∆t δz Kn+1,m 1+1/2+∆t δz Jn+1 o(2.26) •IMC f1=Cn+1,m 1hn+1,m 1+∆t δz Kn+1,m 1+1/2+θn 1−θn+1,m 1+∆t δz Jn+1 o(2.27) Upper boundary aN=−∆t δz2Kn+1,m N−1/2(2.28) bN=Cn+1,m N+∆t δz2Kn+1,m N−1/2(2.29) •IP fN=Cn+1,m Nhn N−∆t δz Kn+1,m N−1/2−∆t δz Jn+1 ∗(2.30) •IMC fN=Cn+1,m Nhn+1,m N−∆t δz Kn+1,m N−1/2+θn N−θn+1,m N−∆t δz Jn+1 ∗(2.31) 24 Chapter 2. Numerical model Lower boundary: Gravity flow/Semi-infinite stratum In order to allow for free gravity flow, the pressure gradient must be zero at the boundary, and because of this, the hydraulic conductivity gradient is also zero. This means h1=h2, hence K1=K2at all times. Finally, flow will be equal to hydraulic conductivity. This can be written in terms of equation (2.37), for i= 2. However, in this particular case it is possible to rewrite coefficients aand bbecause h1=h2, so that h1is not actually solved independently (hence a∗ 2= 0, to remove it from the matrix solution), but later assigned as equal to the solution of h2. a∗ 2hn+1,m+1 1+b∗ 2hn+1,m+1 2+c2hn+1,m+1 3=f2(2.32) a∗ 2= 0 (2.33) b∗ 2=a2+b2(2.34) Impervious stratum In the case that a boundary is an impervious stratum or barrier, the flux at such boundary J1/2or JN+1/2must be imposed as zero. The corresponding coefficients for the first row of the matrix equation are the same as those of the known flow case, with Jo(t) = 0 for the lower boundary or J∗(t) = 0 for the upper boundary. 2.3.4 Convergence and under-relaxation Because Picard iterations are performed in each time step to approximate Cand K, an appropriate convergence criterion is required. The standard is to stop iterations when convergence error hn+1,m+1 i−hn iis less than a specified convergence tolerance ǫ. Huang et al. [22] showed that the standard criterion, although effective, is not particularly efficient in terms of computational time. Huang et al. proposed a θ-based criterion which is computationally more efficient, i.e. θn+1,m+1 i−θn+1,m i. This type of criterion has been also used succesfully by Vanderborght et al. [38], although they found that the θ-based criterion may lead to innacurate matric potential profiles if the value of the C(h) function is small, in which case the standard criterion can be used. Similar experiences are reported by van Dam and Feddes [36]. Both criteria were implemented in this model, although for validation and comparison purposes the h-based criterion has been favored because of the aforementioned observations that the θ-based method may generate errors near saturation. Phoon et al. [28] studied the effects of under-relaxation to accelerate convergence, and conducted a comparison of two under-relaxation methods: UR1, Kn+1,m i=K hn+1,m i+hn i 2!(2.35) and UR2 Kn+1,m i=K hn+1,m i+hn,m−1 i 2!(2.36) Phoon et al. found that, although UR1 is faster than UR2 (and also faster than UR0, i.e, when no under-relaxation is used), it can generate inaccurate results. UR2 improves 2.4. Scheme properties 25 convergence rates compared to no under-relaxation, to a lower degree than UR1 but is much more accurate. In the cases reported in this work, no significant differences in CPU time were found with either methods. 2.4 Scheme properties 2.4.1 Solution method The explicit nature of the EMC and EP schemes imply that no iterations are necessary and each cell can be solved independently in time n+ 1, without the need of a matrix equation, despite the three point finite difference stencil and can be solved by directly evaluating terms in time n. The implicit methods however require the solution of a matrix equation. It is possible to write the schemes (as has been noted already in their description) in a coefficient fashion: aihn+1,m+1 i−1+bihn+1,m+1 i+cihn+1,m+1 i+1 =fi(2.37) It is important to note that term fincludes the water content at time n, and all other terms correspond to time n+ 1, at the miteration. The unknown vector is the matric potential in the entire domain in time n+ 1 and m+ 1 iteration. By writing equation (2.37) from i= 1 to i=N, a linear system can be obtained, with Nequations that can be written in matrix form          b1c10··· 0 a2b2c20. . . 0.........0 . . . 0 aN−1bN−1cN−1 0··· 0aNbN                  hn+1,m+1 1 hn+1,m+1 2. . . hn+1,m+1 N−1 hn+1,m+1 N         =       f1 f2 . . . fN−1 fN        (2.38) This system is clearly tridiagonal and can be solved by Thomas Algorithm, which is an efficient solution of the system. Nevertheless, the implicit schemes require iteration. Hence, for a given time step, the Thomas Algorithm will need to be perfomed as many times as the iterative process requires to converge. In other words, every iteration requires solution of the entire system. 2.4.2 Transition from unsaturated to saturated In neither the EMC nor EP formulations can continuity into the saturated region be achieved. These schemes can solve for the first unsaturated cell iwhich becomes saturated in n+ 1. However, for the next time step the time derivative vanishes in the EMC scheme which in turn vanishes the unknown θn+1 iwhich was to be solved for. In the EP scheme Cn i→0 (as the soil becomes saturated) and the equations are undefined. Hence, calculations must be halted whenever such conditions arise in a time step. Conceptually, the inability of the EMC scheme to solve in saturated conditions is due to the fact that the water-content function is piece-wise defined. It is a smooth function for h≤0, 26 Chapter 2. Numerical model but it is a constant with value θsfor h≥0. Hence, there is no difference in water content for an infinite range of positive pressures. The inability of the EP scheme is of the same nature, although it is through C. In this case, because the state variable to solve for is pressure, from such perspective there is no limitation to variably saturated solutions. But, because the water content function is constant for saturated conditions, then C= 0 in such conditions which undefines the equations. It is in fact the piece-wise definition of θwhich leads to this. Even so, for a one-dimensional case, the solution of the entire profile could be achieved even when saturation conditions arise. Because the suction profile can be determined by the model, the position of the water table can be obtained. Pressure distribution below the water table corresponds to a hydrostatic distribution, thus, the entire pressure profile can be described. Nevertheless, this poses a restriction for the domain of the problem which is indeed solved by the scheme, i.e., the unsaturated region, forcing to define h∈]−∞,0[. This demands that the spatial domain changes size as to coincide with the previous restriction. This is only a computational nuisance, which requires appropriate treatment when coded. The IP and IMC schemes can solve the entire domain even in variably saturated conditions. The issues that impede this in the EMC and EP schemes are not present in the implicit schemes, firstly, because the state variable to solve for is pressure, not water content. Pressure is a continuous function: h∈]−∞,∞[ while the water content function is not. Still, water capacity Cis zero whenever saturation occurs. However, because of the implicit approach, Cis no longer a denominator, but a summand in the non-zero coefficient bi(main diagonal of the coefficient matrix) and in the constant term fi. In summary, to achieve variably saturated solutions, it is necessary to solve for pressure, either directly by discretizing the pressure form (as in the IP scheme) or the mixed form (as in the IMC scheme). However, it is not sufficient to solve for pressure, as can be seen with the EP scheme. To achieve continuity from unsaturated to saturated regimes, the implicit schemes are necessary. 2.4.3 Mass conservation The explicit EMC scheme is a conservative scheme. It solves directly for the conserved state variable. The EP scheme, however, shows very poor conservation properties, despite de fact that it is an explicit scheme. The reason for this is that it solves for pressure head which is not a conserved variable, and relates it to mass conservation through the non-linear function C, which is poorly discretized in time in this scheme, given that the time derivative is evaluated without considering the non-linear relations between C,θand h. Comparison between equations (2.19) and (2.13) shows that only the θterms are present in one and not in the other, and that in equation (2.19) the constant term is dependent on hn+1,m iwhile in equation (2.13) it is dependent on hn i. Note that the difference between the IMC scheme and the IP scheme occurs in the constant term fi. Furthermore, the formulation of boundary conditions changes exactly in the same manner, by applying the same variations in fi. These are the terms responsible for adequate mass conservation [13], whilst all other terms remain identical. The issues that affect the IP scheme which are solved by the IMC scheme, are of the same nature of those responsible for poor conservation in the EP scheme, and are related to the treatment of Cand the time derivative of θ(h). The pressure form solves for a non-conservative state variable, while the water content form, and the mixed, form solve for a conserved variable. From the work of [13] and onwards, the use of the mixed form as in the IMC scheme 2.4. Scheme properties 27 has been widespread. Nevertheless, other approaches have been taken. Rathfelder and Abriola [29] developed mass conservative methods with the h-form by discretizing the hydraulic capacity function with standard chord slope aproximations, but also found that the mixed form, together with an analytical expression for the hydraulic capacity function was computationally more efficient. Phoon et al. [28] also found that under certain under-relaxation techniques, the h-form can be solved with good mass conservation, even in coarse grids. Several authors [13] [22] [28] also emphasize the fact that mass conservation is a necessary condition for an accurate solution, but does not guarantee it. Other factors have severe influence over the accuracy, especially the shape of the hydraulic conductivity function and the effects of discretization on this function. 2.4.4 Stability Stability analysis was performed for the EMC and the IMC schemes only, since basic properties such as mass conservation and accuracy are not well achieved with the EP or the IP scheme. Von Neumann analysis of small perturbations is used. In the case of the EMC scheme, water content perturbations ˜ θof amplitude baround a base value aare done by ˜ θ=a+beiψz (2.39) In the case of the IMC scheme, pressure perturbations ˜ hare studied. ˜ h=a+bei(ψz−ωt)(2.40) Furthermore, because the equation is highly non linear, further assumptions are necessary: Assumption 2.4.1. Water content is a linear function of pressure θ(h) = θo+C(h−ho) (2.41) Assumption 2.4.2. Hydraulic conductivity is a linear function of water content (EMC scheme) or pressure (IMC scheme). K(θ) = Ko+K∗(θ−θo) (2.42) K(h) = Ko+K∗h(h−ho) (2.43) Assumption 2.4.3. The smallest discretized wave which can be observed in a mesh of size δz is of wavelength λ= 4δz. By analyzing the EMC and IMC schemes with these assumptions a linearized analysis can be done, and from such analysis it is possible to obtain some insights of the non-linear stability properties. The reasoning and manipulation can be seen in all extent in Appendix B. The stability condition for the EMC scheme is found to be ∆t≤2δz πν∗1 1 + ǫ∗(2.44) where ν∗=Ko Cwith viscosity units [L2/T] and the dimensionless number ǫ∗=3K∗b Ko . Equation (2.44) shows the charactestic form of a stability condition of a diffusion equation [1] with different coefficients, but nevertheless proportional to the square of mesh resolution, and inversely proportional to a viscosity coefficient. It is important to note that for ǫ∗≪1 is less restrictive than equation (2.44). Note that the slope of the conductivity function, K∗, 28 Chapter 2. Numerical model participates only in ǫ∗. Hence, it is only when ǫ∗is in the order of magnitude of 1 that the K∗affects stability. High values of K∗occur near saturation, specially if using the MG model. Hence, the MMG model is better suited to ensure stability than the MG model. This is consistent with the work by Schaap and van Genuchten [33] and Vogel et al. [39]. The viscocity coefficient ν∗is dependent on C. This allows a simple conclusion. Whenever very dry conditions, or saturation conditions arise, C→0, which implies ν∗→ ∞ ⇒ ∆t→0. In order to further understand the stability properties of EMC to soil parameters, consider the Brooks-Corey model [9] and the Brutsaert equation for conductivity [11]. Then, stability can be written as ∆t=−2δz2 π ω(θs−θr)2 Kshb(θo−θr)θs−θr θo−θr3 2ω    1 1 + 3b2 + 5 2ω(θs−θr)2 θo−θr    (2.45) where the hbis the bubbling pressure (similar to hsin MMG) and ωa fitting parameter that represents pore-size distribution (similar to ˆηin MG and MMG). From (2.45) it can be concluded that for a particular soil, ∆tis inversely proportional to saturation. The higher ωthe more sensible ∆tis to saturation. Conversely, for a particular water content, there is a minimum value of ∆tfor a particular ω. The more saturated the soil is, the least sensitive ∆tis to ω. Saturated conditions result in a ∆twhich varies with ωvery little around the minimum value of ∆t. Hence, saturated, fine-textured soils are very restrictive on time step selection. The analysis of the IMC leads to the conclusion that it is unconditionally stable. This must be considered in context, remembering the linearized analysis from where it derives. Additionally, because the stability analysis was approximated by using Kn, and not Kn+1, a significant part of the non-linearity of the problem might have been lost, specially since an arithmetic mean to obtain Ki±1results in Ki±1=Ko. 2.4.5 Efficiency A priori analysis of the schemes might lead to the conclusion that explicit schemes require less CPU-time since no iterations are necessary. However, because of stability constrains, if the admissible time step is small, CPU time can be much greater than the required CPU time for the implicit schemes. Since efficiency depends on stability constrains, then, the same factors which might lead to small time steps for the EMC scheme affect negatively on EMC efficiency. IMC efficiency depends on the selected time step, but also on soil parameters and water content states. CPU-time requirements for the IMC are dependent on the number of iterations that need to be performed in each time step, times the number of time steps. Larger time steps require more iterations, the question is how many more iterations. An optimal time step might be sought to maximize efficiency. Because iterations intend to linearize K(h), large gradients of hand very non-linear conductivity functions will require more iterations, and hence efficiency is reduced. Nevertheless, it is only through simulation that it can be quantified. Furthermore, time step is likely to be selected according to the desired accuracy and compromising some diffusivity effects, hence, such optimal time step is not investigated in this work. Another factor in efficiency is the efficiency of the algebraic solution of the matrix equation. Because the schemes in this work are only 1D, high efficiency is achieved because of the Thomas Algorithm. However, in more dimensions, the matrix equation is not 2.5. Computation of intercell conductivity Ki±1/229 tridiagonal and efficiency is likely to be severly affected. For such cases, appropriate selection of the alegebraic solver is essential. 2.5 Computation of intercell conductivity Ki±1/2 The issue of selecting an appropriate method to compute the intercell hydraulic conductivity (interblock conductivity, intergrid conductivity, or internode conductivity) has been extensively discussed in the literature and has been identified as a matter of great importance [8], as it can not only affect the quality of the results, but the stability of the numerical model [10]. Several schemes to compute the intercell conductivity have been proposed, analyzed and compared. To illustrate the importance of the method for computing Ki±1/2, consider a discrete domain with constant and small δz. Whenever ∂h ∂z is small between two cells iand i+ 1, the choice of an estimation method for Ki+1/2should not be problematic, as the value will be tightly bounded by Kiand Ki+1 which should be quite similar because of the small difference in pressure (this, however, has been shown not to be true in all cases [3]). Nevertheless, as the gradient of hbecomes larger between two cells, the estimation method becomes important. An inappropriate method, together with the non-linearity of K(h) can lead to large missestimation of Kat the cell interface, thus errors in flow occur. This effect is magnified near saturation, as ∂K ∂h → ∞. This is especially important when solving a problem with boundary conditions that can generate large gradients near the boundaries, because of extreme fixed matric potentials. As fine grids become impractical for large scale or even catchment scale problems, the intercell conductivity estimation method needs to be robust enough to work with relatively coarse grids. The problem of computing intercell conductivity further extends to the estimation of interlayer conductivity in heterogeneous soils [10] [14] [31] and saturated-unsaturated interfaces [27]. Perhaps the most basic method, is the arithmetic mean, Ki±1/2=Ki+Ki±1 2(2.46) The geometric mean has also been proposed, Ki±1/2=pKiKi±1(2.47) The harmonic mean, Ki±1/2=2 1 Ki +1 Ki±1 (2.48) The upstream mean Ki±1/2=         Max(Ki, Ki±1) if ∂h ∂z ≥0 Min(Ki, Ki±1) if ∂h ∂z <0 (2.49) Haverkamp and Vauclin [20] studied several methods and concluded that the geometric mean performs better than other methods. Hornung and Messing [21] showed that the geometric 30 Chapter 2. Numerical model mean performs better than the arithmetic mean. Zaidel and Russo [42] studied the Kirchhoff scheme as well as weighted methods relying on the asymptotic behavior of the conductivity function which for particular cases were reduced to a geometric mean. Van Dam and Feddes [36] used an arithmetic mean although it tends to overestimate infiltration rates (geometric means tend to underestimate it) but concluded that for fine grids the errors generated by the arithmetic mean are smaller than those produced by neglecting hysteresis and spatial soil variability. Gast´o et al. [17] proposed a weighted averages method and found that the arithmetic mean overestimates and the geometric mean underestimates conductivity. Srivastava and Guzman [35] found that integrated conductivity (analitically or by Gaussian integration) provided good results, and also confirmed that the geometric mean outshines the arithmetic and harmonic means, but sometimes even the Gaussian integration method. Another interesting observation is that the upstream conductivity scheme and the harmonic mean scheme provide upper and lower boundaries of the exact solution. Belfort and Lehmann [8] performed simulations with several methods, validating the preference of the geometric mean over the arithmetic and harmonic mean (in particular for large δz), and also finding that for finite elements the geometric mean provides good results and efficieny, but for finite differences, weighted averages can prove better. They concluded that for large nodal spacing arithmetic and upstream means overestimate the wetting front and harmonic and downstream means underestimate it. Vanderborght et al. [38] conducted a set of benchmarking test cases among several codes, concluding that those which use the arithmetic mean predict more dispersed wetting fronts than those obtained with codes that use the geometric mean. Only one code obtained more dispersed fronts with the upstream mean. Warrick [40] who also noted the arithmetic mean and even the geometric mean to be poor estimations, as the geometric mean in some cases greatly underestimated flow. Warrick proposed a weighting scheme which was found to be more accurate but also greatly increased computation time. More recently Baker [2] [3] further analyzed the validity of different means, evaluating if they satisfied mathematical principles (min-max conditions for elliptical value problems) and Darcian flow, and proposed a Darcian mean, i.e., a weighted mean obtained from the spatial distribution of hwhich guarantees darcian flows, finding better accuracy than the geometric mean, but also noted the large computational overhead it requires. Baker showed that by comparing an analytical form of the Darcian mean using the Brooks-Corey model [9], this mean could be reduced to arithmetic, harmonic and geometric means depending on soil parameters, concluding that the arithmetic mean is representative of a nonphysical porous medium, the harmonic relates to an unlikely medium and the geometric mean to clays or rock matrices. Baker’s results show that only the upstream mean did not violate mathematical principles, whilst the arithmetic, harmonic and geometric means showed a great number of violations. Furthermore, Baker showed that traditional means caused non-physical results, which did not occur with the upstream or Darcian mean, and that the main difference between the latter is that the Darcian mean produces sharper wetting fronts and higher peak flows than the upstream mean when space discretization errors occur. Baker’s recommendation is that the upstream mean, because of computer efficiency, is in many cases preferable over the CPU-time-consuming Darcian mean. Because of its simplicity and widespread use, despite the aforementioned studies, the arithmetic mean is included in the model. Nevertheless, following Bakers’ recommendation [3], the upstream mean has also been included, as shown in equation (2.49). A simple way to understand the effects of choosing the arithmetic or upstream mean is to consider the soil shown in Figure 1.3. Consider a downward saturation process, for a cell interface between a saturated cell i= 2 (K→Ks) and a dry boundary cell i= 1 (K→0) with a known and imposed pressure, estimation of the intercell conductivity with an arithmetic mean will result in K1+1/2→Ks/2. If the estimation of the intercell conductivity is done by using the 3.2. Test cases 37 3.2 Test cases A series of test cases were simulated, considering a soil stratum of 100 cm in depth with the following parameters for the MG and MMG models: Ks= 0.00922 cm/s,θs= 0.368 m3/m3, θr= 0.102 m3/m3,α= 0.0335 cm−1, ˆη= 2 (these parameters generate Kand θcurves as shown in figure 1.3). For MMG, hs=−4cm. For all cases a uniform fine grid of ∆t= 1 s, δz = 1 cm was kept as a standard for comparison, and convergence criteria ǫ= 10−7, unless otherwise noted. This is considered to be very strict, in comparison to [22]. For every test case, four simulations were performed to observe sensitivity to particular methods. All test cases are downward processes. A summary of the test cases is presented in table 3.2, and a summary of the setup of the different simulations for each case is presented in table 3.3. Note that only the IMC and EMC schemes were used for these test cases, since the IP and EP schemes were proven inaccurate in validation tests in section 3.1. The EMC scheme was used in those cases in which saturation conditions need not be computed, since the scheme is ineffective when saturation conditions arise. Although cases 6 and 7 appear to be well suited for EMC, when tested, the EMC scheme either was not capable of advancing the drying front (because of initial saturation conditions) or became unstable. Table 3.2: Test Cases Test Case Description Scheme UBC LBC IC 1 Impervious boundaries IMC J∗= 0 cm/s Jo= 0 cm/s h =−20 cm 2 Saturation in semi-infinite soil IMC h= 0 cm ∂h ∂z = 0 h=−100 cm 3 Saturation with water table IMC h= 0 cm cm h = 0 cm h =−100 cm 4 Partial saturation IMC-EMC h=−20 cm h =−50 cm h =−50 cm 5 Full saturation IMC h= 0 cm h =−50 cm h =−50 cm 6 Drying with water table IMC h=−100 cm h = 0 cm h = 0 cm 7 Drying with semi-infinite soil IMC J∗= 0 cm/s ∂h ∂z = 0 h= 0 cm UBC: Upper boundary condition; LBC: Lower boundary condition; IC: Initial condition Table 3.3: Simulation setup Simulation case Mean Soil Model A Arithmetic Mualem-van Genuchten B Upstream Mualem-van Genuchten C Arithmetic Modified Mualem-van Genuchten D Upstream Modified Mualem-van Genuchten 3.2.1 Test Case 1: Impervious boundaries Consider a soil column overlying an impervious stratum (Jo= 0 cm/s) and no infiltration from the surface (J∗= 0 cm/s). Hence, the domain has zero net flow and there is no change in total mass within the soil column. With an initial state of a partially saturated column h(z, t = 0) = −20 cm, the only changes should be the redistribution of water content because of gravity. The simulation was 2 hours long, with results shown every 5 minutes. Simulation results are as expected. No flow enters or exits the domain, and perfect mass balance was obtained in all four simulations. The matric potential profile changes in time, as gravity forces water down, saturating the lower parts of the column, and drying the upper parts until equilibrium (zero flow) matric potential profile is obtained. Theoretically, from 38 Chapter 3. Validation and test cases 0 20 40 60 80 100 −60 −50 −40 −30 −20 −10 0 10 20 30 40 z (cm) h (cm) t = 0 mint = 120 min t = 120 min Case A Case B (a) Simulations A and B 0 20 40 60 80 100 −60 −50 −40 −30 −20 −10 0 10 20 30 40 z (cm) h (cm) t = 0 mint = 120 min t = 120 min Case A Case C (b) Simulations A and C 0 20 40 60 80 100 −60 −50 −40 −30 −20 −10 0 10 20 30 40 z (cm) h (cm) t = 0 mint = 120 min t = 120 min Case B Case D (c) Simulations B and D 0 20 40 60 80 100 −60 −50 −40 −30 −20 −10 0 10 20 30 40 z (cm) h (cm) t = 0 mint = 120 min t = 120 min Case C Case D (d) Simulations C and D Figure 3.5: Results for Test Case 1 equation (1.40), it is clear that this equilibrium requires ∂h ∂z =−1 so that gravitational potential is counteracted. This is verified in the simulation and is easily observed as the uniform slope in the final matric potential profile. Although it is difficult to observe in figure 3.5(a) because of the scale, closer examination shows that in the lower parts of the column (wetting front) in case B, for the same time and depth, shows less hydrostatic pressure than case A. In other words, case B produces a slower wetting front than case A. Conversely, for the drying front in the upper part of the soil column, case B produces slower drying than case A. Nevertheless, as figure 3.5(a) shows, the difference is minimal. Similar behavior is found in 3.5(d). This shows that the upstream mean produces faster wetting fronts compared to the arithmetic mean. Comparison between MG and MMG models (figures 3.5(b) and 3.5(c)) shows clearly that the MMG model generates a faster wetting front, which is due to the fact that maximum conductivity is achieved at lower water contents. 3.2.2 Test Case 2: Downward saturation in semi-infinite soil Initial and boundary conditions were imagined so that the complete saturation process could be observed, allowing gravitational flow in the lower boundary which can be interpreted as having a semi-infinite stratum of the same soil. Initial conditions were h(x, t = 0) = −100 cm. Boundary conditions were ∂h ∂z |x=0,t = 0 and h(x= 100, t) = 0 cm. The simulation was 1 hour long, with results shown every 3 minutes. 3.2. Test cases 39 0 20 40 60 80 100 −100 −80 −60 −40 −20 0 z (cm) h (cm) Case A Case B (a) Simulations A and B 0 20 40 60 80 100 −100 −80 −60 −40 −20 0 z (cm) h (cm) Case A Case C (b) Simulations A and C 0 20 40 60 80 100 −100 −80 −60 −40 −20 0 z (cm) h (cm) Case B Case D (c) Simulations B and D 0 20 40 60 80 100 −100 −80 −60 −40 −20 0 z (cm) h (cm) Case C Case D (d) Simulations C and D Figure 3.6: Results for Test Case 2 Note that although the entire column becomes saturated, pressure head never becomes positive. This is due to the free outflow at the lower boundary by setting the pressure gradient as zero. The effects of the conductivity mean are consistent with those seen in the validation test, showing faster fronts with upstream means. The soil constitutive model generates larger differences in the advancement of the front than mean selection, MMG generates faster fronts than MG. 3.2.3 Test Case 3: Downward saturation with water table This case models a complete saturation process with presence of a fixed water table at the lower boundary. Initial conditions were h(z, t = 0) = −100 cm. Boundary conditions were h(z= 0, t) = 0 cm and h(z= 100, t) = 0 cm. The simulation was 1 hour long, with results shown every 4 minutes. Simulation for Case A became unstable after 20 seconds, and thus no results are shown for clarity, and because of the same reason no results are shown for case C. This case has positive and negative pressure head gradients, which generate a particular scenario for the estimation method for Ki±1/2. The complications are evident in the overshooting effects near the water table boundary, which explain why Simulation A and C failed. Cases B and D do not show overshooting effects, although the negative gradient near the bottom still exists, thus, the arithmetic mean is responsible for the overshooting effects. An additional simulation was performed considering an automatic selection of the mean, in such a way that computation is performed with the arithmetic mean if there is no change in the sign of the 40 Chapter 3. Validation and test cases 0 20 40 60 80 100 −100 −80 −60 −40 −20 0 z (cm) h (cm) Case A Case B (a) Simulations A* and B 0 20 40 60 80 100 −100 −80 −60 −40 −20 0 z (cm) h (cm) Case A Case C (b) Simulations A* and C* 0 20 40 60 80 100 −100 −80 −60 −40 −20 0 z (cm) h (cm) Case B Case D (c) Simulations B and D 0 20 40 60 80 100 −100 −80 −60 −40 −20 0 z (cm) h (cm) Case C Case D (d) Simulations C* and D Figure 3.7: Results for Test Case 3 pressure gradient in the three-point finite difference stencil, and when there is a change in the sign, computation is performed with the upstream mean. Results are reported as case A* and C*. Because cases A and C cannot be properly computed, it is not possible to compare the advance of the wetting front with cases A* and C*. However, it is interesting to note that case A* exhibits a slower wetting front than C* as it is expected. The use of the arithmetic mean as a primary method in cases A* and C* (and the upstream as a correction against overshooting) produces a slightly slower wetting front than cases B and D which are “fully” upstream. 3.2.4 Test Case 4: Downward partial saturation process This case models a saturation process which leads to a partially saturated stationary flow. Initial conditions were h(x, t = 0) = −50 cm. Boundary conditions were h(x= 0, t) = −50 cm and h(x= 100, t) = −20 cm. The simulation time was 3 hours long. Results are shown every 10 minutes. Note that stationary flow occurs around 110 minutes. For this case tests were performes also with the EMC scheme. Results shown are those of simulations with ∆t= 0.7s. The overshooting effect generated by different means can be seen quite clearly. When the arithmetic mean is used, the intercell conductivity in the lower boundary is greatly affected by the unvarying and very low conductivity of the boundary, which in turn is overcompensated by accumulating mass (and hence pressure) in the column. When the upstream mean is 3.2. Test cases 41 0 20 40 60 80 100 −60 −55 −50 −45 −40 −35 −30 −25 −20 −15 −10 z (cm) h (cm) t = 180 min t = 0 min Case A Case B (a) Simulations A and B 0 20 40 60 80 100 −60 −55 −50 −45 −40 −35 −30 −25 −20 −15 −10 z (cm) h (cm) t = 180 min t = 0 min Case A Case C (b) Simulations A and C 0 20 40 60 80 100 −60 −55 −50 −45 −40 −35 −30 −25 −20 −15 −10 z (cm) h (cm) t = 180 min t = 0 min Case B Case D (c) Simulations B and D 0 20 40 60 80 100 −60 −55 −50 −45 −40 −35 −30 −25 −20 −15 −10 z (cm) h (cm) t = 180 min t = 0 min Case C Case D (d) Simulations C and D Figure 3.8: Results for Test Case 4 with IMC used, the lower boundary does not interact with the soil, hence, no overshooting occurs. Nevertheless, in spite of the overshooting effects, the speed of the wetting front is much more sensitive to changes in the consitutive model: MMG produces faster fronts than MG. The small effect produced by selecting upstream or arithmetic means in the speed of the wetting front can be seen in figures 3.8(a) and 3.8(d), which show that upstream produces faster fronts. Figure (3.2.4) shows results for this case using the EMC scheme. The comparisons are the same as those with the IMC scheme, except for figure 3.9(a) which shows results of case A with EMC and IMC. Note that the solution is the same except in the boundary. EMC does not experience overshooting issues. The response of EMC to the use of arithmetic or upstream means is the same as for IMC: upstream means generate faster fronts. The use of MMG generates faster fronts than MG also. In terms of stability, ∆t > 0.7 generated instabilities with the MMG model. With ∆t≈1 instabilities appearead with the MG model. Evaluating equation (2.44) with Ko(h=−20) = 0.0022067, C(h=−20) = 0.103, K∗= 0.06028, with perturbation b= 0.08538 results in ∆tmax = 3.7. Hence, consider the stability number SN =∆treal ∆tmax =0.7 3.7= 0.189 42 Chapter 3. Validation and test cases 0 20 40 60 80 100 −60 −55 −50 −45 −40 −35 −30 −25 −20 −15 −10 z (cm) h (cm) t = 180 min t = 0 min EMC IMC (a) Simulation A vs IMC 0 20 40 60 80 100 −60 −50 −40 −30 −20 −10 z (cm) h (cm) t = 180 min t = 0 min Case A Case B (b) Simulations A and B 0 20 40 60 80 100 −60 −50 −40 −30 −20 −10 z (cm) h (cm) t = 180 min t = 0 min Case A Case C (c) Simulations A and C 0 20 40 60 80 100 −60 −55 −50 −45 −40 −35 −30 −25 −20 −15 −10 z (cm) h (cm) t = 180 min t = 0 min Case B Case D (d) Simulations B and D 0 20 40 60 80 100 −60 −50 −40 −30 −20 −10 z (cm) h (cm) t = 180 min t = 0 min Case C Case D (e) Simulations C and D Figure 3.9: Results for Test Case 4 with EMC 3.2. Test cases 43 3.2.5 Test Case 5: Downward full saturation process This case models the complete saturation process with an imposed pressure head in the lower boundary. Intial conditions were h(x, t = 0) = −50 cm. Boundary conditions were h(x= 0, t) = 0 cm and h(x= 100, t) = −50 cm. Simulation time was 2400 seconds (40 minutes) long. Results are shown every 240 seconds (4 minutes). Simulation B showed convergence issues near saturation, which required to set ǫ= 10−6. Note that stationary flow occured around 24 minutes into the simulation. 0 20 40 60 80 100 −60 −40 −20 0 20 40 60 z (cm) h (cm) t = 40 min t = 0 min Case A Case B (a) Simulations A and B 0 20 40 60 80 100 −60 −40 −20 0 20 40 60 z (cm) h (cm) t = 40 min t = 0 min Case A Case C (b) Simulations A and C 0 20 40 60 80 100 −60 −40 −20 0 20 40 60 z (cm) h (cm) t = 40 min t = 0 min Case B Case D (c) Simulations B and D 0 20 40 60 80 100 −60 −40 −20 0 20 40 60 z (cm) h (cm) t = 40 min t = 0 min Case C Case D (d) Simulations C and D Figure 3.10: Results for Test Case 5 Results for this case are an extension of those of Case 4 into a saturated regime, because the upper boundary condition is set to induce saturation in the entire profile. Conclusions are quite similar to those of Case 4, but in this case it should be noted that the overshooting effect generates a blockage response of the lower boundary, resulting in a fictitious impervious boundary. This shows very clearly that, setting a Dirichlet condition in the lower boundary, and by doing so imposing a conductivity, can dramatically change the pressure profile. It should be noted that it was possible to fully simulate cases A and C, contrary to Case 3 and there was no need to use the automatic selection for averaging K. The effects of the MMG and MG model are as in previous cases, and have no incidence on overshooting which continues to be a conductivity averaging issue. 3.2.6 Test Case 6: Downward drying process with water table This case models a drying process while maintaining a fixed water table in the lower boundary. Initial conditions were set as h(x, t = 0) = 0 cm. Boundary conditions were h(x= 0, t) = 0 cm 44 Chapter 3. Validation and test cases and h(x= 100, t) = −100 cm. This simulation used ∆t= 10sand δz = 2cm. Because simulations A and C became unstable around t= 15 minutes, results are presented every 100 seconds. Simulations B and D were 120 hours long with results shown every 6 hours. Note that stationary flow is reached by the end of the simulation. Comparisons AB and CD provide no interesting information because of the difference in time scales, hence, they are omitted. 0 20 40 60 80 100 −100 −80 −60 −40 −20 0 z (cm) h (cm) t = 15 min t = 0 min Case A Case C (a) Simulations A and C 0 20 40 60 80 100 −100 −80 −60 −40 −20 0 z (cm) h (cm) t = 120 h t = 0 min Case B Case D (b) Simulations B and D Figure 3.11: Results for Test Case 6 Figure 3.11(a) clearly shows the effect of the Dirichlet boundary condition together with the arithmetic mean. It seems clear that the “overdryness” of cell 2 in simulations A and C is the cause of the instability. Because simulations B and D do not show this behavior, it is an indication that the arithmetic mean does not handle well the Dirichlet-type boundary, resulting in overshooting, and eventually, instability. In can be seen from the results that in this drying process, although it is the same soil as in the wetting cases, the use of MMG instead of MG generates faster fronts, but with very little difference in contrast to those generated by MG. 3.2.7 Test Case 7: Downward drying process in semi-infinite soil This test case consists of an initially saturated soil column h(z, t = 0)=0 cm of a semiinfinite soil. Boundary conditions were set as Neumann conditions: no flow from the surface J∗= 0 cm/s and a semi-infinite free draining stratum, which is represented by ∂h ∂z |x=0,t = 0. The simulation was 1 day long, with results shown every hour. Because only Neumann boundary conditions are used, no overshooting issues arise when using the arithmetic mean. Furthermore, the differences generated by selecting the arithmetic mean or the upstream mean are negligible as shown in figures 3.12(a) and 3.12(d). There are far more significant differences when using the MMG or MG constitutive models as can be seen in figures 3.12(b) and 3.12(c), which show that the MMG model generates faster drying fronts. This is curious, since test case 6 resulted in very little differences with respect to both constitutive models. 3.2. Test cases 45 0 20 40 60 80 100 −120 −100 −80 −60 −40 −20 0 z (cm) h (cm) t = 24 h t = 0 h A−MG U−MG (a) Simulations A and B 0 20 40 60 80 100 −120 −100 −80 −60 −40 −20 0 z (cm) h (cm) t = 24 h t = 0 h Case A Case C (b) Simulations A and C 0 20 40 60 80 100 −120 −100 −80 −60 −40 −20 0 z (cm) h (cm) t = 24 h t = 0 h Case B Case D (c) Simulations B and D 0 20 40 60 80 100 −120 −100 −80 −60 −40 −20 0 z (cm) h (cm) t = 24 h t = 0 h Case C Case D (d) Simulations C and D Figure 3.12: Results for Test Case 7 Chapter 4 Conclusions and further research 4.1 Conclusions 1. The EP (Explicit, Pressure form of Richards’ equation) and IP (Implicit, Pressure form of Richards equation) schemes do not converge to correct solutions. They both show very poor conservation properties. 2. The EMC (Explicit, Mixed form of Richards’ equation) is mass conserative and accurate, although incapable of solving saturated conditions, and conditionally stable. 3. The IMC (Implicit, Mixed form of Richards’ equation) scheme, formulated in a similar manner to a pressure form, in such a way as to solve pressure and not water content is appropriate to solve variably saturated flow, with adequate mass conservation, and reasonable efficiency. 4. EMC and IMC scheme approximate correctly the solution to Richards’ equation. Differences in the solutions were only observed near Dirichlet boundaries, where de EMC scheme generates smooth transitions, while de IMC scheme results in discontinuities or even overshooting. 5. The EMC scheme can be computationally more efficient than IMC in certain cases (particularly for the same ∆t). However, in other cases EMC can be severely less efficient because of stability constrains. 6. The stability analysis of the IMC scheme and the test cases show that the scheme is inconditially stable. 7. The stability analysis of the EMC scheme shows that it is conditionally stable, in a way similar to that of a diffusion equation. Valditation and tests cases confirm such dependence, although, because of the non-linearities of the equation, the expression for maximum time step is not precise, but a guideline of stability requirements. 8. Richards’ equation in the water content form cannot solve variably saturated flow problems, only unsaturated flow problems. 9. Richards’ equation in the mixed form cannot solve -per sevariably saturated flow problems. Appropriate discretization of the time derivative in order to obtain a discretized mixed form which solves for pressure is necessary. A.3. Implicit Presure based scheme 53 Cn i hn+1 i−hn i ∆t=1 δz Kn i+1/2 hn i+1 −hn i δz −Kn i−1/2 hn i−hn i−1 δz +Kn i+1/2−Kn i−1/2 hn+1 i=∆t Cn iδz2Kn i+1/2hn i+1 −∆t Cn iδz2Kn i+1/2+Kn i−1/2hn i+∆t Cn iδz2Kn i−1/2hn i−1 +∆t Cn iδz Kn i+1/2−Kn i−1/2+hn i(A.2) A.3 Implicit Presure based scheme From Richards’ equation in pressure form, with a backward Euler scheme for time derivatives and a centered finite-difference scheme for spatial derivatives, Cn+1 hn+1 −hn ∆t−∂ ∂z "K(hn+1)∂hn+1 ∂z + 1#= 0 In order to formulate the Picard iteration, let δm=hn+1,m+1 −hn+1,m where m+ 1 is the computed iteration and mthe previous iteration. Cn+1,m δm ∆t+Cn+1,m hn+1,m −hn ∆t−∂ ∂z Kn+1,m ∂δm ∂z =∂Kn+1,m ∂z +∂ ∂z Kn+1,m ∂hn+1,m ∂z  Spatial discretization with finite differences, using ∂δm ∂z i±1/2≈δm i−δmi±1 δz and evaluating hydraulic conductivity in between cells yields Cn+1,m i δm i ∆t+Cn+1,m i hn+1,m i−hn i ∆t−1 δz2hKn+1,m i+1/2δm i+1 −δm i−Kn+1,m i−1/2δm i−δm i−1i= 1 δz Kn+1,m i+1/2−Kn+1,m i−1/2+1 δz2hKn+1,m i+1/2hn+1,m i+1 −hn+1,m i−Kn+1,m i−1/2hn+1,m i−hn+1,m i−1i Expanding all terms, Cn+1,m i δm i ∆t+Cn+1,m i hn+1,m i−hn i ∆t−1 δz2Kn+1,m i+1/2δm i+1+1 δz2Kn+1,m i+1/2δm i+1 δz2Kn+1,m i−1/2δm i−1 δz2Kn+1,m i−1/2δm i−1= 1 δz Kn+1,m i+1/2−Kn+1,m i−1/2+1 δz2Kn+1,m i+1/2hn+1,m i+1 −1 δz2Kn+1,m i+1/2hn+1,m i−1 δz2Kn+1,m i−1/2hn+1,m i+1 δz2Kn+1,m i−1/2hn+1,m i−1 Grouping terms by spatial index i−1, i, and i+ 1 −1 δz2Kn+1,m i−1/2δm i−1−1 δz2Kn+1,m i−1/2hn+1,m i−1+Cn+1,m iδm i ∆t+Cn+1,m i hn+1,m i ∆t−Cn+1,m i hn i ∆t +1 δz2Kn+1,m i+1/2δm i+1 δz2Kn+1,m i−1/2δm i+1 δz2Kn+1,m i+1/2hn+1,m i+1 δz2Kn+1,m i−1/2hn+1,m i −1 δz2Kn+1,m i+1/2δm i+1 −1 δz2Kn+1,m i+1/2hn+1,m i+1 = 1 δz Kn+1,m i+1/2−Kn+1,m i−1/2 54 Appendix A. Numerical schemes formulations Expanding δmwith its definition, −1 δz2Kn+1,m i−1/2hn+1,m+1 i−1+1 δz2Kn+1,m i−1/2hn+1,m i−1−1 δz2Kn+1,m i−1/2hn+1,m i−1 +Cn+1,m ihn+1,m+1 i ∆t−Cn+1,m ihn+1,m i ∆t+Cn+1,m i hn+1,m i ∆t−Cn+1,m i hn i ∆t +1 δz2Kn+1,m i+1/2hn+1,m+1 i−1 δz2Kn+1,m i+1/2hn+1,m i +1 δz2Kn+1,m i−1/2hn+1,m+1 i−1 δz2Kn+1,m i−1/2hn+1,m i+1 δz2Kn+1,m i+1/2hn+1,m i+1 δz2Kn+1,m i−1/2hn+1,m i −1 δz2Kn+1,m i+1/2hn+1,m+1 i+1 +1 δz2Kn+1,m i+1/2hn+1,m i+1 −1 δz2Kn+1,m i+1/2hn+1,m i+1 = 1 δz Kn+1,m i+1/2−Kn+1,m i−1/2 Grouping and rearranging, results in −∆t δz2Kn+1,m i−1/2hn+1,m+1 i−1+Cn+1,m i+∆t δz2Kn+1,m i+1/2+Kn+1,m i−1/2hn+1,m+1 i+−∆t δz2Kn+1,m i+1/2hn+1,m+1 i+1 = ∆t δz Kn+1,m i+1/2−Kn+1,m i−1/2+Cn+1,m ihn i A.4 Implicit Mixed Conservative scheme From Richards’ equation in the mixed form, with a backward Euler scheme for time derivatives and a centered finite-difference scheme for spatial derivatives, θn+1 −θn ∆t−∂ ∂z "K(hn+1)∂hn+1 ∂z + 1#= 0 In order to formulate the Picard iteration, let δm=hn+1,m+1−hn+1,m where m+1 is the computed iteration and mthe previous iteration. Invoking Taylor’s polynomial to approximate θn+1,m+1 in terms of the derivative C=∂θ ∂h and δm, θn+1,m +Cn+1,mδm−θn ∆t−∂ ∂z Kn+1,m ∂δm ∂z =∂Kn+1,m ∂z +∂ ∂z Kn+1,m ∂hn+1,m ∂z  Spatial discretization with finite differences, using ∂δm ∂z i±1/2≈δm i−δmi±1 δz and evaluating hydraulic conductivity in between cells yields θn+1,m i+Cn+1,m iδm i−θn i ∆t−1 δz2hKn+1,m i+1/2δm i+1 −δm i−Kn+1,m i−1/2δm i−δm i−1i= 1 δz Kn+1,m i+1/2−Kn+1,m i−1/2+1 δz2hKn+1,m i+1/2hn+1,m i+1 −hn+1,m i−Kn+1,m i−1/2hn+1,m i−hn+1,m i−1i Expanding all terms, θn+1,m i−θn i ∆t+Cn+1,m iδm i ∆t−1 δz2Kn+1,m i+1/2δm i+1+1 δz2Kn+1,m i+1/2δm i+1 δz2Kn+1,m i−1/2δm i−1 δz2Kn+1,m i−1/2δm i−1= 1 δz Kn+1,m i+1/2−Kn+1,m i−1/2+1 δz2Kn+1,m i+1/2hn+1,m i+1 −1 δz2Kn+1,m i+1/2hn+1,m i−1 δz2Kn+1,m i−1/2hn+1,m i+1 δz2Kn+1,m i−1/2hn+1,m i−1 A.4. Implicit Mixed Conservative scheme 55 Grouping terms by spatial index i−1, i, and i+ 1 −1 δz2Kn+1,m i−1/2δm i−1−1 δz2Kn+1,m i−1/2hn+1,m i−1 +Cn+1,m iδm i ∆t+1 δz2Kn+1,m i+1/2δm i+1 δz2Kn+1,m i−1/2δm i+1 δz2Kn+1,m i+1/2hn+1,m i+1 δz2Kn+1,m i−1/2hn+1,m i −1 δz2Kn+1,m i+1/2δm i+1 −1 δz2Kn+1,m i+1/2hn+1,m i+1 = 1 δz Kn+1,m i+1/2−Kn+1,m i−1/2−θn+1,m i−θn i ∆t Expanding δmwith its definition, −1 δz2Kn+1,m i−1/2hn+1,m+1 i−1+1 δz2Kn+1,m i−1/2hn+1,m i−1−1 δz2Kn+1,m i−1/2hn+1,m i−1 +Cn+1,m ihn+1,m+1 i ∆t−Cn+1,m ihn+1,m i ∆t+1 δz2Kn+1,m i+1/2hn+1,m+1 i−1 δz2Kn+1,m i+1/2hn+1,m i +1 δz2Kn+1,m i−1/2hn+1,m+1 i−1 δz2Kn+1,m i−1/2hn+1,m i+1 δz2Kn+1,m i+1/2hn+1,m i+1 δz2Kn+1,m i−1/2hn+1,m i −1 δz2Kn+1,m i+1/2hn+1,m+1 i+1 +1 δz2Kn+1,m i+1/2hn+1,m i+1 −1 δz2Kn+1,m i+1/2hn+1,m i+1 = 1 δz Kn+1,m i+1/2−Kn+1,m i−1/2−θn+1,m i−θn i ∆t Grouping and rearranging, results in −∆t δz2Kn+1,m i−1/2hn+1,m+1 i−1+Cn+1,m i+∆t δz2Kn+1,m i+1/2+Kn+1,m i−1/2hn+1,m+1 i+−∆t δz2Kn+1,m i+1/2hn+1,m+1 i+1 = Cn+1,m ihn+1,m i+∆t δz Kn+1,m i+1/2−Kn+1,m i−1/2+θn i−θn+1,m i Appendix B Stability Analysis B.1 EMC Scheme Recall equation (1.40), the mixed form of Richards’ equation: ∂θ(h) ∂t =∂ ∂z "K(h)∂h ∂z + 1#=−∂J ∂z In order to perform the stability analysis consider the following: •Conductivity is assumed as linear function of water content: K(θ) = Ko+K∗(θ−θo) •Water content is assumed as linear function of water head: θ(h) = θo+C(h−ho). Consider a soil column with an initially uniform moisture distribution in the entire depth. Let ˜ θbe a perturbation of the initial moisture distribution, described by ˜ θ=a+beiψz (B.1) In consequence, ∂˜ θ ∂z =∂ ∂z a+beiψz=ibψeiψz Because θ=θ(h) ∂˜ θ ∂h ∂h ∂z =C∂h ∂z =ibψeiψz Rearranging, ∂h ∂z =ibψ Ceiψz (B.2) From the definition of the flux J, and considering equation (B.2) J=−K∂h ∂z + 1=−Kibψ Ceiψz −K Introducing the linear form of K, J=−(Ko+K∗θ−K∗θo)ibψ Ceiψz + 1 B.1. EMC Scheme 57 Hence, flux ˜ Jfor the perturbation θ=˜ θis ˜ J=−(Ko−K∗θo)ibψ Ceiψz + 1−K∗a+beiψzibψ Ceiψz + 1 Regrouping, ˜ J=−Ko+K∗θo−K∗a−(Ko−K∗θo+K∗a+K∗b)ibψ Ceiψz−K∗ ib2ψ Ce2iψz Because θo=afor perturbation conditions, ˜ J=−Ko−(Ko+K∗b)ibψ Ceiψz |{z } ˜ J1 −K∗ ib2ψ Ce2iψz | {z } ˜ J2 (B.3) A first order Taylor expansion of ˜ J(ξ) in the neigbourhood of ξ=ϕis ˜ J(ξ) = ˜ J(ϕ) + ∂˜ J(ϕ) ∂ξ (ξ−ϕ) + O(ξ2) Consider ξ=eiψz and ϕ=ξ(z= 0) = 1. The derivative required for the Taylor expansion is ∂˜ J ∂ξ =∂˜ J1(ξ) ∂ξ +∂˜ J2(ξ) ∂ξ where ∂˜ J1(ξ) ∂ξ =−(Ko+K∗b)ibψ C and ∂˜ J2(ξ) ∂ξ =−2K∗ ib2ψ Cξ Hence, ∂˜ J ∂ξ =−(Ko+K∗b)ibψ C−2K∗ ib2ψ Cξ In consequence, the Taylor polynomial, neglecting O(ξ2), is evaluated as ˜ J=−Ko−(Ko+K∗b)ibψ C−2K∗ ib2ψ C−h(Ko+K∗b)ibψ C+ 2K∗ ib2ψ Cieiψz −1 Finally, the perturabated flux is ˜ J=−Ko−(Ko+ 3K∗b)ibψ Ceiψz (B.4) Consider the EMC scheme, written in terms of flux J, θn+1 i=θn i−∆t δz Jn i+1/2−Jn i−1/2(B.5) 58 Appendix B. Stability Analysis Consider that in time nthe perturbation ˜ θand the perturbation flow ˜ Joccur, hence θn=˜ θ and Jn=˜ J. Hence, Jn i±1/2=˜ J(z±δz). Consequently, by substituting equations (B.1) and (B.4) into equation (B.5) yields θn+1 i=a+beiψzi+∆t δz "Ko+ (Ko+ 3K∗b)ibψ Ceiψ(z+δzi/2)# −∆t δz "Ko+ (Ko+ 3K∗b)ibψ Ceiψ(z−δzi/2)# Cancelling terms and rearranging, θn+1 i=a+beiψzi+∆t δz (Ko+ 3K∗b)ibψ Ceiψzeiψδz/2−∆t δz (Ko+ 3K∗b)ibψ Ceiψze−iψδz/2 Grouping, θn+1 i=a+beiψzi1 + ∆t δz (Ko+ 3K∗b)iψ Ceiψδz/2−e−iψδz/2(B.6) Equation (B.6) is the expression which describes how a perturbation in an an initially uniform water content profile progresses in time. In order for the scheme to be stable, the perturbation cannot amplify itself. Hence,  1 + ∆t δz (Ko+ 3K∗b)iψ Ceiψδz/2−e−iψδz/2 2 ≤1 Because of the relation of the complex exponential function to tringonometric functions,  1 + ∆t δz (Ko+ 3K∗b)iψ C2isin ψδz 2 1 ≤1 Hence,  1−2ψ C ∆t δz (Ko+ 3K∗b) sin ψδz 2 2 ≤1 Consider that the smallest wave length that can be observed in a uniform space discretization of size δz is λ= 4δz, which implies that the the wave number satisifies ψ=π 2δz , which implies sin ψδ 2= sin π 4=√2 2. Hence,  1−π∆t Cδz2Ko+ 3K∗b 2 ≤1 Expanding the squared modulus, 1 + π2∆t2 C2δz4Ko+ 3K∗b2−2π∆t Cδz2Ko+ 3K∗b≤1 Solving for ∆t, ∆t≤2 π δz2C Ko+ 3K∗b B.1. EMC Scheme 59 Consider the definition of a dimensionless number ǫ∗, ǫ∗=3K∗b Ko (B.7) And consider the following defintion of parameter ν∗with viscosity units, ν∗=Ko C(B.8) Then, the stability criterion is ∆t≤2 π δz2 ν∗1 1 + ǫ∗(B.9) Consider equation (B.9), where ǫ∗≪1, then ∆t≤2 π δz2 ν∗ Note that ν∗has viscosity units L2 Tand that the stability criterion resembles the well-known stability criterion for an explicit centered finite difference scheme for the 1D diffusion equation ∆t≤δx2 2αwhere αis a constant diffusion coefficient [1]. However, ν∗is not a constant coefficient, and the stability criterion depends on Koand C. Effects of ν∗and ǫ∗ To further investigate stability, consider the Brooks-Corey model [9]: θ=hb hω (θs−θr) + θr(B.10) where hbis the bubbling pressure and ωa fitting parameter that represents pore-size distribution. Hence, its derivative is C=∂θ ∂h =−ωhω b(θs−θr)h−(ω+1) (B.11) Consider the conductivity function as K=Ksθ−θr θs−θr2+ 5 2ω as suggested by Brutsaert [11]. Then, K∗=∂K ∂θ =2 + 5 2ωKs (θs−θr)5 2ω (θ−θr)5 2ω+1 (B.12) Hence, ν∗=Ko C=− Ksθo−θr θs−θr2+ 5 2ω ωhω b(θs−θr)h−(ω+1) Which is ν∗=−Ks ωhb θo−θr (θs−θr)2θo−θr θs−θr3 2ω 60 Appendix B. Stability Analysis On the other hand, ǫ∗=3K∗b Ko = 3b 2 + 5 2ωKs (θs−θr)5 2ω (θ−θr)5 2ω+1 Ksθ−θr θs−θr2+ 5 2ω Rearranging, ǫ∗= 3b2 + 5 2ω(θs−θr)2 θo−θr Finally, the stability criterion can be expressed as ∆t=−2δz2 π ω(θs−θr)2 Kshb(θo−θr)θs−θr θo−θr3 2ω    1 1 + 3b2 + 5 2ω(θs−θr)2 θo−θr     From this expression it can be concluded that, for a particular soil, ∆tis inversely proportional to saturation. The higher ωthe more sensible ∆tis to saturation. Conversely, for a particular humidity, there is a minimum value of ∆tfor a particular ω. The more saturated the soil is, the least sensitive ∆tis to ω. Saturated conditions result in a ∆twhich varies with ωvery little around the minimum value of ∆t. B.2 IMC Scheme Consider the perturbation of an initial state ˜ h=a+bei(ψz−ωt)(B.13) Hence, hn+1 i=a+bei(ψzi−ωtn+1) Assume that the perturbation does not suffer changes in wave length in space, hence: hn+1 i+1 =a+bei(ψ(zi+δz)−ωtn+1 )=a+bei(ψzi−ωtn+1)ei(ψδz)=a−aeiψδz +hn+1 ieiψδz In summary, hn+1 i+1 =a−aeiψδz +hn+1 ieiψδz (B.14) And in a similar way, hn+1 i−1=a−ae−iψδz +hn+1 ie−iψδz (B.15) Let hydraulic conductivity and water content be linear functions of pressure, hence K=Ko+Kh(h−ho) (B.16) θ=θo+C(h−ho) (B.17) B.2. IMC Scheme 61 Consider the IMC scheme as described by equation (2.15) −∆t δz2Kn+1,m i−1/2hn+1,m+1 i−1+Cn+1,m i+∆t δz2Kn+1,m i+1/2+Kn+1,m i−1/2hn+1,m+1 i +−∆t δz2Kn+1,m i+1/2hn+1,m+1 i+1 =Cn+1,m ihn+1,m i+∆t δz Kn+1,m i+1/2−Kn+1,m i−1/2+θn i−θn+1,m i Recall that iterations in this scheme respond mainly to mass-balance issues because of the linearization of C. Hence, for small amplitude perturbations of h(and consequently of Cand K) consider time m+ 1 as n+ 1 for hand time mas time nfor Kand C. Note that for the scheme to converve in a single time step it is necessary that θm≈θm+1. Hence, the scheme may be rewritten as −∆t δz2Kn i−1/2hn+1 i−1+Cn i+∆t δz2Kn i+1/2+Kn i−1/2hn+1 i −∆t δz2Kn i+1/2hn+1 i+1 =Cn ihn i+∆t δz Kn i+1/2−Kn i−1/2+θn i−θn+1 i(B.18) Substituting equations (B.14) and (B.15) in equation (B.18) and the linear definition of θ, −∆t δz2Kn i−1/2a−ae−iψδz +hn+1 ie−iψδz+Cn i+∆t δz2Kn i+1/2+Kn i−1/2hn+1 i −∆t δz2Kn i+1/2a−aeiψδz +hn+1 ieiψδz=Cn ihn i+∆t δz Kn i+1/2−Kn i−1/2 +θo+Cn i(hn i−ho)−θo−Cn i(hn+1 i−ho) Grouping, ∆t δz2hn+1 i2δz2 ∆tCn i−Kn i−1/2e−iψδz +Kn i+1/2+Kn i−1/2−Kn i+1/2eiψδz −a∆t δz2hKn i−1/21−e−iψδz+Kn i+1/21−eiψδzi= 2Cn ihn i+∆t δz Kn i+1/2−Kn i−1/2 (B.19) Because of the linear definition of K, Ki±1/2=Ko+Kh(hi−ho) + Ko+Kh(hi±1−ho) 2=2Ko+Khhi+Khhi±1−2Khho 2 Kn i±1/2=Ko−Khho+Kh 2hn i+hn i+±1 Furthermore, hn i±1=a−ae±iψδz +hn ie±iψδz. Hence, Kn i±1/2=Ko−Khho+Kh 2hn i+a−ae±iψδz +hn ie±iψδz Note that because of the perturbation analysis, a=hoand hn i=ho, which results in Kn i±1/2=Ko Let Cn i=Ch. In consequence, equation (B.19) becomes ∆t δz2Kohn+1 i2δz2 ∆t Ch Ko + 2 −eiψδz +e−iψδz−hoKo ∆t δz2h2−eiψδz +e−iψδzi= 2Chho 62 Appendix B. Stability Analysis By means of the identities between de complex exponential function and trigonometric functions, ∆t δz2Kohn+1 iδz2 ∆t Ch Ko + 1 −cos (ψδz)−hoKo ∆t δz2[1 −cos (ψδz)] = Chho The amplification factor is G=hn+1 i ho , hence Gδz2 ∆t Ch Ko + 1 −cos (ψδz)−1 + cos (ψδz) = Ch δz2 Ko∆t Solving for G G= 1 + Ch δz2 Ko∆t−cos (ψδz) 1 + Ch δz2 Ko∆t−cos (ψδz) = 1 (B.20) This result implies that for small perturbations, the IMC method is unconditionally stable. If the analysis is performed for the first iteration, i.e., θm+1 =θnthe same conclusion is obtained.