# Modeling

The THYRSIS hydrogeological flow and transport modeling is based of the coupling of a 1D vertical flow and transport model in the unsaturated zone and a 2D horizontal flow and transport model in the saturated zone. The saturated 2D modeling is based on a steady state flow field, which is created independently for each site.

This steady state flow field can be imported from an external calculation or created directly in THYRSIS by defining a 2D mesh with the GMSH mesh builder [GR09] and determining the flow characteristics - mainly the permeability field - with the HYNVERSE inverse process [RJ08], [RT16].

## Steady state groundwater flow

Steady state groundwater flow must be defined first, because it is used for transitional flow and transport calculations in saturated zone. It is also needed to define the characteristics of the 1D unsaturated columns, namely :

- the column thicknesses,

- the saturated permeability in the columns.

Steady state groundwater flow is a 2D calculation whose parameters are :

- infiltration, which is constant in steady state, but can be spatially heterogeneous,

- groundwater thickness, variable on the whole modeled domain,

- permeability field.

The output are :

- potential (or piezometric) field, i.e the altitude of the water table in case of unconfined groundwater.

- Darcy’s velocity field.

Futhermore, it is the same mesh that is used for flow and transport in saturated zone.

All these informations are saved in spatialite databases.

Steady state groundwater flow is a complex hydrogeological modeling process. The model is built using a great amount of data concerning the topography and the geology (in order to determine the aquifer vertical extension), the meteorology (in order to estimate the infiltration) and the hydrogeology. These data are the result of bibliographic studies and borehole analysis on the site and its neighborhood, because the domain extent is often wider than the studied site. The permeability field definition is a hard task, mostly resolved using inversion calculation, coupling geostatistics and optimization process (HYNVERSE [RJ08], [RT16]).

## Groundwater flow in transitional regime

Groundwater flow in transitional regime is based on the model built in permanent regime, adding a storage capacity value, equal to the effective porosity in case of unconfined groundwater.

The flow simulation in transitional regime needs the knowledge of the infiltration variation over time. This information is provided as an external file, according to the appropriate syntax (cf. [User files]).

## Flow and transport in unsaturated zone

Unsaturated zone modeling is based on 1D vertical columns creation. There are as many columns as injection zones, and their thicknesses are calculated by substracting the steady state piezometric level to the topographic altitude.

Flow is calculated using the Richards equation, with the Van Genuchten formula [VG80], modified by [IVB06]:

\[\begin{split} S_{e}=\begin{cases} \frac{1}{S_{c}}\left[1+\left(\alpha h\right)^{n}\right]^{-m} & h>h_{e}\\ 1 & h\leq h_{e} \end{cases} \end{split}\]
\[\begin{split} K=\begin{cases} K_{s}\sqrt{S_{e}}\left[\frac{1-\left(1-\left(S_{e}S_{c}\right)^{1/m}\right)^{m}}{1-\left(1-S_{c}^{1/m}\right)^{m}}\right]^{2} & S_{e}<1\\ K_{s} & S_{e}\geq1 \end{cases} \end{split}\]

with :

- \(h\) : pressure potential (L)

- \(h_{e}\) : input pressure (L)

- \(K\) : permeability (\(L.T^{-1}\))

- \(K_{s}\) : saturated permeability (\(L.T^{-1}\))

- \(S_{e}=\frac{\theta-\theta_{r}}{\theta_{s}-\theta_{r}}\) : reduced water content, or effective saturation

- \(S_{c}=\left[1+\left(\alpha h_{e}\right)^{n}\right]^{-m}\)

- \(\theta\) : water content

- \(\theta_{r}\) : residual water content

- \(\theta_{s}\) : saturated water content

- \(m=1-\frac{1}{n}\)

- \(n\) and \(\alpha\) (\(L^{-1}\)) are fitting parameters.

There is a one to one connection between water content \(\theta\) and the saturation permeability \(s\) using \(\omega\) : \(\theta=s\omega\) .

When \(h_{e}=0\) , then \(S_{c}=1\) and we get the original Van Genuchten formula.

### Mesh building of unsaturated columns

The columns are split in centimetric thickness elements, on one hand to keep numeric stability even for low dispersivity value, and on the other hand to garanty accurate mass balance.

More precisely, at the top of the column, the thickness of the first element is :

\[ min\left(1cm,\frac{e_{ZNS}}{100}\right) \]

where \(e_{ZNS}\) is the thickness of the unsaturated column, expressed in centimeters.

The first element’s thickness is mostly equal to 1 cm, except for columns with a thickness less than 1 m. Then element’s thicknesses increase with a geometric progression until the middle of the columns, to finally dicrease in symmetry to reach the bottom of the column. The geometric progression is deduced from the ratio between the thickness of the middle column element and the thickness of the top column element. The ratio is equal to 5 by default and can be configured by the user (coef_maillage in [simulations] table).

Furthermore, the elements thickness is constrained to be less than \(0.1/\alpha\) (van Genuchten coefficient) in order to garanty numeric stability in case of high \(\alpha\) value.

### Unsaturated zone parameters

The unsaturated zone differs from the saturated zone, because the flow and the transport are processed in a single simulation, which starts with a steady state flow calculation, followed by a transitional flow and transport calculation.

The transport model parameters in the unsaturated zone are :

- kinematic porosity \(\omega_{c}\) (fracture and/or matrix pore system, according to the studied site),

- total porosity \(\omega_{t}\),

- volumetric mass density \(\rho_{s}\),

- \(\alpha\) and n van Genuchten coefficients,

- maximal and residual saturations,

- longitudinal dispersivity (transverse dispersivity doesn’t exist in 1D),

- partition coefficient \(K_{d}\) ,

- chemical solubility \(L_{s}\) (cf. [Mass injection]).

Among these parameters, the total porosity and the volumetric mass density are used to calculate the time delay coefficient R, according to the following formula :

\[ R=1+\frac{\left(1-\omega_{t}\right)}{\omega_{c}}\rho_{s}K_{d} \]

As for saturated zone, transitional infiltration simulations in unsaturated zone require the knowledge of the infiltration variation over time, which is provided as an external file, according to the appropriate syntax (cf. appendix [User files]).

### Diffuse injection on the whole domain

It is possible to simulate a diffuse injection on the whole domain, either by simulating as many unsaturated vertical columns as mesh nodes, or by simulating columns groups (composed of similar thickness and permeability) to reduce the amount of simulations. In this case, two parameters must be defined ([simulations] table) :

- pas_echant_zns is the column sampling step. It defines the column size of a same group. If the step is defined to 5 m, groups will be defined as follow [0,5m], [5m,10m], [10m,15m], etc. The more the step is high, the more the number of groups is low.

- rapport_max_permeabilite is the maximal ratio of permeability tolerated in a same group of columns, by default equal to 10.

## Transport in saturated zone

Transport in saturated zone is based on the flow model in [permanent] regime or in [transitional] regime.

A injection in the groundwater can be :

- a flow, a concentration or a mass, defined via the interface, if the migration in unsaturated zone has not to be modeled (cf. [injection]),

- a flow, defined from flows calculated at the bottom of unsaturated columns. In this case, the flow is injected into the meshes (OpenFOAM) or nodes (METIS) related to each injection zones, in order to cover the whole injection surface.

Transport model paramaters in saturated zone are :

- porosity (fracture and/or matrix pore system, according to the studied site),

- longitudinal and transverse dispersivities,

- chemical retention rate.

In case of a flow in transitional regime, the storage capacity coefficient is set to the effective porosity.

## Injection characteristics

Chemical injection is defined by the following characteristics :

- coordinates (x, y) of the injection zone center,

- depth,

- area S concerned by the injection,

- flow F, expressed in kg/s or Bq/s,

- injection duration \(\Delta T\).

For convenience, the injection can be defined as a concentration or mass, as described below. A leak can also be modeled. To define it, a water volume is required.

### Concentration injection

In case of an injection defined by a concentration `C`, flow `F` is calculated from :

\[ F=ISC \]

where `I` is the rain water infiltration, defined with the flow model in [permanent] regime.

### Mass injection

In case of an injection defined by a mass `M`, flow `F` and duration \(\Delta T\) are calculated using the solubility limit \(L_{s}\) of the chemical :

\[ F=ISL_{s} \]
\[ \Delta T=\frac{M}{F} \]

### Leak

It is possible to add a water injection to the chemical injection by defining the water volume injected \(V\), in addition to the flow \(F\) or the concentration \(C\). In that case, injection velocity is calculated from :

\[ v=\frac{V}{S \Delta T} \]

This velocity can’t be greater than the saturated permeability \(K_{s}\). If it is the case, more precisely if :

\[ v > 0.9 K_{s} \]

injection duration is modified by :

\[ \Delta T=\frac{V}{0.9 S K_{s}} \]

and a message is emited to notify that it is advised to increase the injection area by setting :

\[ S=\frac{V}{0.9 K_{s} \Delta T} \]

If concentration \(C\) is provided, then flow is equal to :

\[ F=vSC \]

The following spreadsheet summarizes the flow calculation, according to the available data :

Data | \(V=0\) | \(V>0\)
\(F, \Delta T\) | \(F\) | \(v=min(\frac{V}{S \Delta T}, 0.9 K_{s})\), \(\Delta T = max(\Delta T, \frac{V}{0.9 S K_{s}})\)
\(C, \Delta T\) | \(F=ISC\) | \(v=min(\frac{V}{S \Delta T}, 0.9 K_{s})\), \(\Delta T = max(\Delta T, \frac{V}{0.9 S K_{s}})\), \(F=vSC\)
\(M, L_{s}\) | \(F=ISL_{s}\) | No calculation

## Uncertainty analysis

THYRSIS provides the opportunity to create several flow and transport parameter sets, according to standard probability laws (uniform, log-uniform, normal, log-normal), in order to get confidence intervals on simulated variables or to create probability maps according to a threshold.

The number of simulations must be greater than 1 and a probability law and min/max values must be defined for each parameter. The parameter sets are generated by [Latin Hypercube], and one simulation is realized for each parameter set. The simulations can be parallelized according to the settings in the [Preferences] . Finally, THYRSIS aggregates the results in order to calculate mean and standard deviation for each simulated variable.

Graphics display confidence intervals or probability maps.

## Dual porosity model

Dual-porosity flow model, as explained by Gerke & Van Genuchten [GVG93], is available in THYRSIS, but currently only when using METIS. This model is based on the assumption that the Richards equation applies simultaneously in the matrix pore system and in the fracture pore system. For a one dimensional vertical system, flow equations in a elementary volume are :

In fracture pore system :

\[ C_{f}\frac{\partial h_{f}}{\partial t}=\frac{\partial}{\partial z}\left(K_{f}\frac{\partial h_{f}}{\partial z}-K_{f}\right)-\frac{\Gamma_{w}}{w_{f}}-S_{f} \]

In matrix pore system :

\[ C_{m}\frac{\partial h_{m}}{\partial t}=\frac{\partial}{\partial z}\left(K_{m}\frac{\partial h_{m}}{\partial z}-K_{m}\right)-\frac{\Gamma_{w}}{1-w_{f}}-S_{m} \]

with :

\[ C_{f,m}=S_{w}S_{s}^{f,m}+\omega_{f,m}\frac{\partial S_{w}}{\partial h} \]
\[ \Gamma_{w}=\alpha_{w}\left(h_{f}-h_{m}\right) \]
\[ \alpha_{w}=\alpha_{w}^{*}K_{a} \]
\[ \alpha_{w}^{*}=\frac{\beta}{a^{2}}\gamma_{w} \]

and :

\(K_{f}\) and \(K_{m}\) (\(L.T^{-1}\)) : fracture and matrix pore system permeabilities,

\(h_{f}\) and \(h_{m}\) (L) : fracture and matrix pore system hydraulic head,

\(\Gamma_{w}\) (\(T^{-1}\)) : exchange term describing water transfert between fracture and matrix pore system,

\(S_{f}\) and \(S_{m}\) (\(T^{-1}\)) : source terms in fracture and matrix pore system,

\(w_{f}\) : volume fraction of the fracture pore system in the total system,

\(C_{f}\) and \(C_{m}\) (\(L^{-1}\)) : specific water capabilities in fracture and matrix pore system,

\(S_{w}\) : water saturation rate,

\(S_{s}\) (\(L^{-1}\)) : specific storage capacity coefficient in fracture and matrix pore system,

\(\omega_{f}\) and \(\omega_{m}\) : porosity in fracture and matrix pore system,

\(K_{a}\) (\(L.T^{-1}\)) : transition permeability between fracture and matrix pore system,

\(a\) (L) : mean half-distance between fractures,

\(\beta\) : structure factor, that is the geometry of the fracture system (beta=3 in case of a set of orthogonal and regularly spaced fractures),

\(\gamma_{w}\) : empiric factor equal to 0,4.