## 2 Materials and Methods
### Groundwater Modelling  
In this work, the MODFLOW [47] code was chosen to simulate the groundwater flow of the Tebessa-Morsott alluvial aquifer system. The MODFLOW code allows the calculation of the hydraulic head on each mesh of a three-dimensional domain by the finite difference method [48], using seeking to solve the diffusivity equation with partial derivatives of flows by a combination of the continuity equation and Darcy's law. The representative elementary volume adapted to the flow in the aquifer must consider the flow domain over the entire wetted height between the coasts z1 and z2. Level z1 represents the impermeable substratum of the aquifer, and level z2 represents either the impermeable cover of a confined aquifer or the piezometric surface of an unconfined aquifer (as in our case), where the coast is identified as the head h. The governing equations for 2D or 3D groundwater flow are based on water mass balance and [PERSON]'s law [49]:  
\[\left[\frac{\partial}{\partial x}\left(K_{x}\frac{\partial h}{\partial x} \right)+\frac{\partial}{\partial y}\left(K_{y}\frac{\partial h}{\partial y} \right)+\frac{\partial}{\partial z}\left(K_{z}\frac{\partial h}{\partial z} \right)\right]=\ S_{s}\frac{\partial h}{\partial t}+q_{s} \tag{2}\]  
By integrating and taking into account: \(\partial h/\partial z=0\) ([PERSON]'s hypothesis):  
\[\int_{z1}^{z2}\frac{\partial}{\partial x}(K_{x}\frac{\partial h}{\partial x} )dz+\int_{z1}^{z2}\frac{\partial}{\partial y}\left(K_{y}\frac{\partial h}{ \partial y}\right)dz=\int_{z1}^{z2}S_{s}\frac{\partial h}{\partial t}dz+\int_ {z1}^{z2}qdz \tag{3}\]  
Again, assuming that z1 and z2 vary in relation to x and y:  
\[\frac{\partial}{\partial x}((\int_{z1}^{z2}K_{X}dx)\frac{\partial}{\partial x })+\frac{\partial}{\partial y}((\int_{z1}^{z2}K_{y}dy)\frac{\partial}{ \partial y})=(\int_{z1}^{z2}s_{s}dz)+(\int_{z1}^{z2}qdz) \tag{4}\]  
In this work, the modeling steps were:  
1. Mesh and geometric criteria: the mesh is an important step that should be adapted to the problem. The size of the modeled domain depends on several factors: (i) The desired precision in the calculation; (ii) the number and distance of water points (boreholes and wells); and (iii) the shape of the piezometric surface in relation to the limit of the aquifer (modeled domain).  
2. Boundary conditions: the different types of boundary conditions of the model are inspired by the geology of the region. They make it possible to delimit the domain within which the model will be responsible for calculating the hydraulic head of each mesh from the initial loads introduced as input. It is important to ensure that the boundaries imposed on the model create the same effect as its natural boundaries: (i) The model is limited vertically by a substratum--this is a zero-flow limit; (ii) the lateral limits are also considered  
\begin{table}
\begin{tabular}{c c} \hline \hline
**DRASTIC Vulnerability** & **Vulnerability Index** \\ \hline Low & 23–100 \\ Average & 101–140 \\ High & 141–200 \\ Very high & 200–226 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Evaluation classes of vulnerability degree using the DRASTIC model.  
as zero-flow (impermeable) limits; and (iii) at the surface, the water table is free, and its load is expressed by an infiltration rate--this is a condition of distributed flow over the surface of the aquifer roof.