CEMRACS 2019 Coupling model of underground flow and pollution - PowerPoint PPT Presentation

CEMRACS 2019 Coupling model of underground flow and pollution transport using a Finite volume scheme Ayoub Charhabil charhabil@math.univ-paris13.fr F. Benkhaldoun Paris 13 University August 8th, 2019 A.Charhabil Paris 13 IRSN 1 / 45

Introduction In Nature, It is more complex ! Figure: Porous media A.Charhabil Paris 13 IRSN 2 / 45

Methematical model Mathematical Model The problem is fully descriebed by the following system of equations : Richards Equation : ∂θ ( h , x , y ) = ∇ . [ K ( h , x , y ) ∇ h ] + ∇ . K ( h , x , y ) + Q s ∂ t Transport Equation : ∂θ C + ∇ . ( qC ) = ∇ . ( θ D ∇ C ) ∂ t , A.Charhabil Paris 13 IRSN 3 / 45

Methematical model 3D Richards Equation : ∂θ ( h , x , y , z ) = ∇ . [ K ( h , x , y , z ) ∇ h ] + ∇ . K ( h , x , y , z ) + Q s ∂ t With : h : head water θ : volumetric water content K : hydraulic conductivity Q s :source term A.Charhabil Paris 13 IRSN 4 / 45

Methematical model 3D Transport Equation : � ∂θ C ∂ t + ∇ . ( qC ) = ∇ . ( θ D ∇ C ) θ D = λ | q | + θ D m τ with : C : solute concentration | q | :Darcy velocity λ longitudinal lenth pore/solide θ : volumetric water content A.Charhabil Paris 13 IRSN 5 / 45

Methematical model 1D Richards Equation : The Richards equation takes 3 forms : The θ -Form : ∂θ ( h ) = ∂ � � ∂θ ∂ z + ∂ K �� D ( θ ) ∂ t ∂ z ∂ z The Mixed-Form : ∂θ ( h ) = ∂ � � ∂ h �� K ( θ ) ∂ z − 1 ∂ t ∂ z The h-Form : C ( h ) ∂ h ∂ t = ∂ � � ∂ h �� K ( h ) ∂ z − 1 ∂ z A.Charhabil Paris 13 IRSN 6 / 45

Methematical model 1D Richards Equation : The θ -Form : ∂θ ( h ) = ∂ � � ∂θ ∂ z + ∂ K �� D ( θ ) ∂ t ∂ z ∂ z With : D = K ∂ h [ L 2 T − 1 ] ∂θ Advantages: Conservation form by construction Mass balance is improved significantly Rapid convergence Inconvinients: limited to unsaturated conditions (In saturation D is infinite !) Limited to homogenous soil ( θ can be not continuous across interfaces separating the layers !) A.Charhabil Paris 13 IRSN 7 / 45

Methematical model 1D Richards Equation : The Mixed-Form : � � ∂ h �� ∂θ ( h ) = ∂ K ( θ ) ∂ z − 1 ∂ t ∂ z Advantages: Mass Conservation / Mass balance Applicable to both saturated and unsaturated soil Applicable to heterogenous soil Inconvinients: Acceptable numerical solutions not always guarenteed A.Charhabil Paris 13 IRSN 8 / 45

Methematical model 1D Richards Equation : The h-Form : C ( h ) ∂ h ∂ t = ∂ � � ∂ h �� K ( h ) ∂ z − 1 ∂ z With : C ( h ) = ∂θ (C:capillary capacity) ∂ h Advantages: Applicable to both saturated and unsaturated soil Applicable heterogenous soil Very close to the physical model Less complicated to implement Inconvinients: Poor preservation of mass balance Relatively slow convergence A.Charhabil Paris 13 IRSN 9 / 45

Coupling system 1D-Coupling : We choose the h-Form ! Coupling of h-Form of Richards and Transport equations in 1D : � ∂ h  C ( h ) ∂ h ∂ � �� ∂ t = K ( h ) ∂ z − 1  ∂ z      � ∂ h  q = − ∂ � �� K ( h ) ∂ z − 1 ∂ z      − ∂ qC ∂θ C ∂ � θ D ∂ C �  ∂ t =  ∂ z ∂ z ∂ z K ? C ? θ ? A.Charhabil Paris 13 IRSN 10 / 45

The parameters of the physical model: The Brooks-Corey model: The Hydraulic conductivity is : � K ( h ) = K s [ θ h − θ r θ s − θ r ] 3+2 / n If h < h d K ( h ) = K s If h � h d With : K s Hydraulic conductivity in saturation h d is the bubbling or air entry pressure head (L) and is equal to the pressure head to desaturate the largest pores in the medium n = 1 − 1 / m ,m parameters linked to the soil structure A.Charhabil Paris 13 IRSN 11 / 45

The parameters of the physical model: The Brooks-Corey model: The Capillary capacity is taken as followed : � C ( h ) = n θ s − θ r | h d | ( h d h ) n +1 If h < h d C ( h ) = 0 If h � h d NB : The Capillary capacity is always positive ! θ s water content in saturation θ r residual water content n et m parameters linked to the soil structure A.Charhabil Paris 13 IRSN 12 / 45

The parameters of the physical model: The Brooks-Corey model: the volumetric water content is taken as followed : � θ = θ r + ( θ s − θ r ) h d If h < h d h θ = θ s If h � h d With : θ s water content in saturation θ r residual water content n et m parameters linked to the soil structure A.Charhabil Paris 13 IRSN 13 / 45

The parameters of the physical model: The van Genuchten Model: Capillary capacity is taken as followed : � n ∗ m ∗ a | h | d θ C ( h ) = If h < 0 (1+( a ∗| h | ) n ) 1+ m C ( h ) = 0 If h � 0 NB : The Capillary capacity is always positive ! With: d θ = θ s − θ r S ∗ The specific volumetric storativity a ,n et m parameters linked to the soil structure A.Charhabil Paris 13 IRSN 14 / 45

The parameters of the physical model: The van Genuchten Model: We introduice the saturation S e as followed : 1 S e = (1 + a n | h | n ) m The Hydraulic conductivity is :  √ S e (1 − � 1 e ) m  K ( h ) = K s 1 − S If h < 0 m K ( h ) = K s h � 0 If  With: K s Hydraulic conductivity in saturation a ,n et m parameters linked to the soil structure NB : Hydraulic conductivity is always positive ! A.Charhabil Paris 13 IRSN 15 / 45

The parameters of the physical model: The van Genuchten Model: There is a ”relationship” between θ and S e , and it’s formulated this way : In saturation case : S e = θ − θ r θ s − θ r In non saturation case : S e = 1 A.Charhabil Paris 13 IRSN 16 / 45

Coupling system Finite Volumes Scheme : General form we use the following FV-schemes : Explicite : h n +1 = h n j − r [ φ ( h n j , h n j +1 ) − φ ( h n j , h n j − 1 )] j Implicite : h n +1 j − r [ φ ( h n +1 , h n +1 j +1 ) − φ ( h n +1 , h n +1 = h n j − 1 )] j j j with φ the numerical flux (In our case, we consider 2 numerical flux adequate to our study :flux of ROE or Lax-Frederick ) A.Charhabil Paris 13 IRSN 17 / 45

Coupling system Richards Equation :Numerical scheme The Explicite Finite volume scheme : r h n +1 = h n j )( φ n j +1 / 2 − φ n j + j − 1 / 2 ) j C ( h n with : r = ∆ t ∆ z φ : numerical flux for h K ( h n h n j +1 − h n j +1 / 2 ) j φ n j +1 / 2 = − j +1 / 2 ) ( − 1) θ ( h n ∆ z For the explicit version of our model The stability condition is : ∆ t ≤ CFLInfC ∗ ∆ z 2 2 Maxk A.Charhabil Paris 13 IRSN 18 / 45

Coupling system Richards Equation :Numerical scheme The implicite Finite volume scheme : r h n +1 = h n )( φ n +1 j +1 / 2 − φ n +1 j + j − 1 / 2 ) j C ( h n +1 j with : r = ∆ t ∆ z φ : numerical flux for h K ( h n +1 h n +1 j +1 − h n +1 j +1 / 2 ) φ n +1 j j +1 / 2 = − j +1 / 2 ) ( − 1) θ ( h n +1 ∆ z A.Charhabil Paris 13 IRSN 19 / 45

Coupling system Transport Equation :Numerical scheme We use an upwind scheme (1st order) : 2 C n +1 = C n +1 − r ∗ V ∗ ( fluxS n j − fluxS n ∗ ( DiffS n j − DiffS n j − 1 ) + r ∗ j − 1 ) j j θ n j + θ n j +1 With : qS n j + qS n V = j +1 2 q n qS n j j = θ n j � fluxS n j = C n If q � 0 j fluxS n j = C n q < 0 If j +1 dz ∗| q n j | DiffS n ∗ ( C n j +1 − C n j = j ) / dz Pe A.Charhabil Paris 13 IRSN 20 / 45

Coupling system Numerical resultats Soil parameters for our test case : Figure: The soil A.Charhabil Paris 13 IRSN 21 / 45

Coupling system Numerical resultats Soil parameters for our test case : case Sand Infiltration through homogenous Domaine lenth L=100 cm K s = 0 . 00922 cm / s , θ s = 0 . 368, θ r = 0 . 102, a = 0 . 0335 cm − 1 Parameters case 1 ∆ z = 0 . 2 cm , Pe = 0 . 2 case 2 ∆ z = 2 cm , Pe = 20 case 3 ∆ z = 2 cm , Pe = 200 case 2 Infiltration Through hetergenous soil Domaine lenth L=100 cm K s = 0 . 000151 cm / s , θ s = 0 . 4686, θ r = 0 . 106, a = 0 . 03104 cm − 1 Parameters (Clay) A.Charhabil Paris 13 IRSN 22 / 45

Coupling system Numerical resultats Head water in 5 and 10 days : temps en jours: 5 temps en jours: 10 0 0 -1 -1 -2 -2 -3 -3 Pression effective h Pression effective h -4 -4 -5 -5 -6 -6 -7 -7 -8 -8 -9 -9 -10 -10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 profondeur z profondeur z A.Charhabil Paris 13 IRSN 23 / 45

Coupling system Numerical resultats : Solute concentration Solute concentration in 5 and 10 days : temps en jours: 5 temps en jours: 10 1 1 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 Solute s Solute s 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 profondeur z profondeur z A.Charhabil Paris 13 IRSN 24 / 45

Coupling system Numerical resultats :Hydraulic conductivity Hydraulic conductivity in 5 and 10 days: # 10 -7 temps en jours: 5 # 10 -6 temps en jours: 10 7 1 0.9 6 0.8 5 0.7 0.6 4 K(h) K(h) 0.5 3 0.4 0.3 2 0.2 1 0.1 0 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 profondeur z profondeur z A.Charhabil Paris 13 IRSN 25 / 45