Mass-Conservative Finite Volume Methods for the Richards equation - PowerPoint PPT Presentation

Mass-Conservative Finite Volume Methods for the Richards’ equation Lecturer: G. Manzini 1 1 IMATI - CNR, Pavia Siena 2006 (Joint collaboration with Stefano Ferraris, Univ. of Turin, Italy)

R + , Richards’ equation: mixed θ − ψ formulation R d , computational domain • The solution ψ ( x , t ) satisfies ∂θ ( ψ ) � � − div K ( ψ ) ∇ ( ψ + z ) = s , Ω × in ∂ t • Ω ⊂ • ψ , hydraulic pressure head; • θ ( ψ ) , volumetric water content; • K ( ψ ) , hydraulic conductivity tensor; • z , vertical coordinate (positive upward) • s, production/consumption source term;

R + , Richards’ equation: head formulation • The solution ψ ( x , t ) satisfies η ( ψ ) ∂ψ � � ∂ t − div K ( ψ ) ∇ ( ψ + z ) = s , in Ω × where • η ( ψ ) = ∂θ ( ψ ) /∂ψ is the general storage term. • Remark: the two formulations are theoretically equivalent, but leads to different numerical discretizations.

R + , R + , To complete the model. . . . . . we must assign: • a proper set of boundary conditions : Γ D × ψ = g D , Dirichlet BCs: on Γ N × − n · K ( ψ ) ∇ ( ψ + z ) = g N , Neumann BCs: on where ∂ Ω = Γ D ∪ Γ N , and g D and g N are, respectively, Dirichlet and Neumann boundary data; • an initial solution : ψ ( x , t = 0 ) = ψ 0 ( x ) x ∈ Ω , for • the characteristic relations for θ ( ψ ) , K ( ψ ) , and η ( ψ ) depending on soil nature: Van Genuchten , e.g. HYDRUS-2D, 1999, Brooks-Corey , Haverkamp (1977), and many other models. . .

Extension to multi-layered soils Multi-layered soils are taken into consideration by assuming that Ω be partitioned into a set of soil zones (sub-surface layers) Ω l , l = 1 , . . . , L such that � Ω = Ω l , l = 1 ,..., L showing different constitutive relations: θ ( ψ ) | Ω l = θ l ( ψ ) , K ( ψ ) | Ω l = K l ( ψ ) , η ( ψ ) | Ω l = η l ( ψ ) .

The finite volume framework � � Let T h = be a suitable triangulation of Ω ; T Let us integrate over any control volume T ∈ T h to obtain � � � � � ∂θ ( ψ ) dV − K ( ψ ) ∇ ( ψ + z ) dV = s dV ; div ∂ t T T T then, apply the Gauss-Green divergence theorem � � � ∂θ ( ψ ) � � n · dV − K ( ψ ) ∇ ( ψ + z ) dS = s dV , ∂ t T ∂ T T introduce cell averages and split face contributions to flux balance � � � d � � � n · θ ( ψ ) dV − | e | K ( ψ ) ∇ ( ψ + z ) dS = s dV . dt T ∂ T T e ∈ ∂ T

Finally, let us introduce the cell-average approximation of the pressure field, e.g. ψ h ≡ { ψ T } , and the numerical flux, � � G ij ( ψ h ) = | e ij | n ij · K ij G ⋄ ij ( ψ h ) + ζ , where ζ = ∇ z is the versor along Z , � � K ij ( ψ h ) = 1 1 1 1 K ( ψ j ) + , 2 K ( ψ i ) and G ⋄ ij ( ψ h ) is the face-based gradient on the diamond D ij that is given by ij ( ψ h ) n ij + ψ β − ψ α ij ( ψ h ) = G ( n ) t ij . G ⋄ | e ij |

The normal component G ( n ) ij ( ψ h ) of the face-based gradient G ⋄ ij ( ψ h ) is 1 � G ( n ) w ( n ) ij , k ψ k + g ( n ) ij ( ψ h ) = ij |D ij | k � � w ( n ) where the weights are obtained by the Least Squares ij , k reconstruction algorithm. Rewriting � � � | e ij | n ij · K ij G ⋄ G ( ψ ) ψ + g ( ψ ) i = ij ( ψ h ) , e ij ∈ ∂ T i � � � | e ij | n ij · K ij ζ , h ( ψ ) i = e ij ∈ ∂ T i we finally obtain the vector formulation d θ ( ψ ) � � � � i − G ( ψ ) ψ + g ( ψ ) + h ( ψ ) i = s i ( ψ ) . � dt

Time discretization • Crucial issue : we are advancing in time θ ( ψ ) , but we need ψ for flux balancing: d θ ( ψ ) � � � � + G ( ψ ) ψ + g ( ψ ) + h ( ψ ) i = s i ( ψ ) . � dt i � �� � � �� � � �� � non-linear gravity term non-linear flux term function of ψ • the simplest approach is to switch to a discrete form of the head-based formulation by using d θ i ( ψ i ) ≈ η ( ψ i ) d ψ i dt . dt

• Using the implicit Euler F.D. scheme for the time derivative give us ) ψ n + 1 − ψ n η ( ψ n + 1 i i i ∆ t � � G ( ψ n + 1 ) ψ n + 1 + g ( ψ n + 1 ) + h ( ψ n + 1 ) i = s ( ψ n + 1 ) . − • This non-linear algebraic problem can be solved by using a fixed-point iterative scheme (Picard): � � � � ψ n + 1 , m + 1 = η ( ψ n + 1 , m − ∆ t G ( ψ n + 1 , m ) ) diag i � � η ( ψ n + 1 , m ) ψ n + ∆ t s i ( ψ n + 1 , m ) + g ( ψ n + 1 , m ) + h ( ψ n + 1 , m )

• The (so-called) delta-form is obtained by solving for δ m + 1 , m � � � ψ n + 1 , m + 1 − ψ n + 1 , m � � � η ( ψ n + 1 , m ) − ∆ t G ( ψ n + 1 , m ) = diag � �� � = δ m + 1 , m � � s ( ψ n + 1 , m ) + g ( ψ n + 1 , m ) + h ( ψ n + 1 , m ) +∆ t � �� ψ n + 1 , m − ψ n � η ( ψ n + 1 , m ) + ∆ t G ( ψ n + 1 , m ) ψ n + 1 , m − diag • Iterate on m to get ψ n + 1 , m + 1 from ψ n + 1 , m starting from ψ n + 1 , 0 = ψ n until � � δ m + 1 , m � � ≤ TOL and, finally, set ψ n + 1 , m + 1 = ψ n + 1 , m + δ m + 1 , m This approach is easy, but shows very poor mass conservation!

What is a mass-conservative discretization? Following M. Celia et al. , W.R.R. 1990, let us define:       total additional total mass in total mass in  =  −  , mass in the the domain at the domain at    domain any time t > 0 initial time t = 0 � flux balance integrated from initial � � � total net flux = into the domain time t = 0 to the current time t � MASS BALANCE RATIO � = total additional mass in the domain total net flux into the domain Clearly, it must be MASS BALANCE RATIO =1.

How does the head-based approach perform? Let us consider, for example, the vertical column infiltration problems proposed by M. Celia et al. W.R.R. 1990: Case 1 Data: [bottom] 0 ≤ ζ ≤ 40 [top] column range: − 61 . 5 cm bottom pressure: − 20 . 7 cm top pressure: − 61 . 5 cm for 0 ≤ ζ < 40 initial pressure: from T = 0 to T = 360 s. running time: Constitutive relationships (Haverkamp et al.): θ ( ψ ) = α ( θ s − θ r ) α + | ψ | β + θ r A K ( ψ ) = K s A + | ψ | γ ( α, β, γ, θ s , θ r , K s and A are known parameters.)

Column infiltration problem: Case 1 Mass balance ratio for different time steps ∆ t 1 Head-based 2D FV scheme (M. & Ferraris 2004) Head-based 1D FD scheme (Celia et.al., 1990) 0.95 Mass Balance Ratio 0.9 0.85 0.8 0.75 0 100 200 300 400 Time-step size (sec)

How does the head-based approach perform? Case 2 Data: [bottom] 0 ≤ ζ ≤ 100 [top] column range: − 1000 cm bottom pressure: − 120 cm top pressure: initial pressure: − 1000 cm for 0 ≤ ζ < 100 running time: from T = 0 to T = 1 day. Constitutive relationships (Van Genuchten et al.): θ s − θ r θ ( ψ ) = 1 + ( ε | ψ | n ) m + θ r � � − m � 2 1 − ( ε | ψ | ) n − 1 � 1 + ( ε | ψ | n ) K ( ψ ) = K s � m � 1 + ( ε | ψ | n ) 2 ( m , n , ε, θ s , θ r , and K s are known parameters.)

Column infiltration problem: Case 2 Mass balance ratio for different time steps ∆ t 1 Head-based 2D FV scheme (M. & Ferraris, 2004) Head-based 1D FD scheme (Celia et.al. 1990) 0.9 Mass Balance Ratio 0.8 0.7 0.6 0.5 0.4 0 1000 2000 3000 4000 Time-step size (sec)

How can we improve this behavior? • Main Idea: M. Celia et al. observed that mass conservation may be lost in the head-based formulation because of a poor approximation of the time derivative of θ ( ψ ) : � � � � ∂θ ( ψ ) n + 1 1 ∂θ ( ψ ) n + 1 dV ≈ d n + 1 � � � ≈ dt θ ( ψ i ( t )) � � � ∂ t | T i | ∂ t i i T i i ) ψ n + 1 − ψ n ≈ η ( ψ n i i + O ( | δψ | ) ∆ t • Better strategy: develop θ ( ψ ) to second-order in ψ -variation in the cell-average integral � 1 θ n + 1 , m + 1 θ ( ψ n + 1 , m + 1 = ) dV i i | T i | T i

How can we improve this behavior? A straightforward calculation gives � 1 θ n + 1 , m + 1 θ ( ψ n + 1 , m + 1 = ) dV i i | T i | T i � 1 � � ψ n + 1 , m + ( ψ n + 1 , m + 1 − ψ n + 1 , m = θ ) dV i i i | T i | T i � �� � = δ m + 1 , m i � � � � ∂θ � n + 1 , m � 1 � � θ ( ψ n + 1 , m ψ n + 1 , m + 1 − ψ n + 1 , m + O ( | ψ | 2 ) � � = ) + dV � i i i | T i | ∂ψ i T i ≈ θ n + 1 , m + η ( ψ n + 1 , m ) δ m + 1 , m + O ( | δψ | 2 ) . i i

• Let us start again from � � � d θ � i − G ( ψ ) ψ + g ( ψ ) + h ( ψ ) i = s i ( ψ ) , � dt • and discretize directly the time derivative of θ ( ψ ) to get θ n + 1 − θ n � � G ( ψ n + 1 ) ψ n + 1 + g ( ψ n + 1 ) + h ( ψ n + 1 ) i i = s i ( ψ n + 1 ) . i − ∆ t • Use Picard iterations for solving the non-linear problem θ n + 1 , m + 1 − θ n i i ∆ t � � G ( ψ n + 1 , m ) ψ n + 1 , m + 1 + g ( ψ n + 1 , m ) + h ( ψ n + 1 , m ) i = s i ( ψ n + 1 , m ) . −