## 2 Materials and Methods
### Description of the Numerical Model  
#### 2.2.1 Code  
SUFT3D [39; 40] is the finite element numerical code we used to develop the groundwater numerical model. This code solves the groundwater flow equation (Equation (1)) based on a mixed formulation of [PERSON]'s equation proposed by [PERSON] et al. [41] using the control volume finite element (CVFE):  
\[\frac{\partial\theta}{\partial t}=\underline{\
abla}.\underline{K}(\theta) \underline{\
abla}h+\underline{\
abla}.\underline{K}(\theta)\underline{ \
abla}z+q, \tag{1}\]  
where \(t\) is the time [T], \(\theta\) is the water content [-], \(z\) is the elevation [L], \(h\) is the pressure head [L], \(q\) is a source/sink term [T\({}^{-1}\)], and \(\underline{K}\) is the hydraulic conductivity tensor [LT\({}^{-1}\)] defined as  
\[\underline{K}=K_{r}\underline{K}_{S}, \tag{2}\]  
where \(\underline{K}_{S}\) is the saturated permeability tensor [LT\({}^{-1}\)], and \(K_{r}\) is the relative hydraulic conductivity [-] that varies from a value of 1 for full saturation to a value of 0 when the water phase is considered immobilized [42]. In the partially saturated zone, the value of \(K_{r}\) evolves according with the following equations [43]:  
\[\theta=\theta_{r}\frac{(\theta_{s}-\theta_{r})}{(h_{b}-h_{a})}(h-h_{a}), \tag{3}\]  
\[K_{r}(\theta)=\frac{\theta-\theta_{r}}{\theta_{s}-\theta_{r}} \tag{4}\]  
where \(\theta_{r}\) is the residual water content [-], \(\theta_{s}\) is the saturated water content [-], \(h_{a}\) is the pressure head at which the water content is just lower than the saturated one [L], and \(h_{b}\) is the pressure head at which the water content is the same as the residual one [L]. The \(K_{r}\) varies linearly between the unsaturated and saturated zones as can be observed in Equations (3) and (4). This adopted linearity for defining the transition between saturated and unsaturated zones does not alter the results of the model, because this work is focused on processes that occurred only in the saturated portion of the soil, while this contributed to mitigate the convergence errors that are common when non-linear expressions are used. The main reason we choose SUFT3D is because it has certain capabilities specifically designed for modelling underground mines, improving the realism and the results of the groundwater numerical model. Specifically, underground cavities were simulated as linear reservoirs using the hybrid finite element mixing cell (HFEMC) method [39; 40] implemented in the SUFT3D code [44; 45; 46]. This method combines physically-based and spatially distributed models as well as black-box models. The domain can be divided into different subdomains depending on their characteristics with specific behaviors assigned depending on their nature. Groundwater processes through unmined areas, which were modeled with finite elements, were computed according to the flow equation in variably saturated porous media (Equation (1). Single mixing cells were used to discretize the underground cavities (i.e., chambers) that are modelled as linear reservoirs, which is similar to the box model approaches that consider a mean hydraulic head for the whole cell. The groundwater exchange between the domains modelled as linear reservoirs and those modelled as porous medium varies linearly [47] and is governed by the following internal dynamic Fourier boundary condition (BC) [40]:  
\[Q_{i}=\alpha^{\prime}A\left(h_{aq}-h_{ur}\right), \tag{5}\]  
where \(h_{aq}\) is the piezometric head in the aquifer [L], \(h_{ur}\) is the hydraulic head in the underground reservoir [L], \(Q_{i}\) is the exchanged flow [L\({}^{3}\)T\({}^{-1}\)], \(A\) is the exchange area [L\({}^{2}\)], and \(\alpha^{\prime}\) is the exchange coefficient [T\({}^{-1}\)]. It is important to highlight that the water velocity inside the mixing cell is not considered; however, the influence of this particularity on the results was minimized by using different linear reservoirs for modelling the different chambers.