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 Richard's equation proposed by Celia et al. [41] using the control volume finite element (CVFE):

$$\frac{\partial \theta}{\partial t} = \underline{\nabla} \cdot \underline{\underline{K}}(\theta) \underline{\nabla} h + \underline{\nabla} \cdot \underline{\underline{K}}(\theta) \underline{\nabla} z + q\_{\prime} \tag{1}$$

where *t* is the time [T], *θ* is the water content [-], *z* is the elevation [L], *h* is the pressure head [L], *q* is a source/sink term [T−<sup>1</sup> ], and *K* is the hydraulic conductivity tensor [LT−<sup>1</sup> ] defined as

$$
\underline{\underline{K}} = \mathcal{K}\_r \underline{\underline{K}}\_{S'} \tag{2}
$$

where *K<sup>S</sup>* is the saturated permeability tensor [LT−<sup>1</sup> ], and *K<sup>r</sup>* 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<sup>r</sup>* 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 *θ<sup>r</sup>* is the residual water content [-], *θ<sup>s</sup>* is the saturated water content [-], *h<sup>a</sup>* is the pressure head at which the water content is just lower than the saturated one [L], and *h<sup>b</sup>* is the pressure head at which the water content is the same as the residual one [L]. The *K<sup>r</sup>* 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–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' A \left( h\_{aq} - h\_{ur} \right), \tag{5}$$

where *haq* is the piezometric head in the aquifer [L], *hur* is the hydraulic head in the underground reservoir [L], *Q<sup>i</sup>* is the exchanged flow [L3T −1 ], *A* is the exchange area [L<sup>2</sup> ], and *α* 0 is the exchange coefficient [T−<sup>1</sup> ]. 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.

SUFT3D was also used because it allows adopting virtual connections that are essential for the purpose of this investigation. Virtual connections, that are also named "by-pass", allow the establishment of hydraulic connections between non-adjacent subdomains that are modeled as linear reservoirs. Virtual connections are defined by a first-order transfer equation (Equation (5) that can be switched off or on according to the hydraulic head difference between the two connected subdomains:

$$Q\_{\upsilon r} = \mathfrak{a}\_{\upsilon r} \left( h\_{SDj} - h\_{SDi} \right) \tag{6}$$

where *hSDj* and *hSDi* are the hydraulic heads inside each of the connected linear reservoirs, *Qvr* is the flow between reservoirs, and *αvr* is the exchange coefficient of the virtual connection [L2T −1 ]. Virtual connections allow constraining the maximum and minimum hydraulic heads into the underground reservoir. They are crucial for the development of the model, as it is not possible to anticipate the water exchanges and, therefore, when the underground reservoir will be full or empty.

Consequently, it is not possible to predict if a pumping or discharge can be carried out. Two types of virtual connections are implemented. The first one extracts immediately and automatically the discharged water when the underground reservoir is full. The virtual connection is switched off for most of the time, and it is only switched-on when the underground reservoir is completely full. At this moment, a value of 10<sup>6</sup> *m*2/*d* is adopted for *αvr* to force the surplus of discharged water to be extracted. The second virtual connection is used to avoid the hydraulic head being lower than the bottom of each chamber. In this case, the connection is switched on most of the time (*αvr* = 10<sup>6</sup> *m*2/*d* to allow the hydraulic connection between chambers). Only when the hydraulic head is at a lower elevation than the bottom of a chamber is the virtual connection deactivated, thus, disconnecting individually each chamber from the rest of the underground reservoir if the hydraulic head is too low.
