2.2.2. Characteristics of the Model

The numerical model is a squared domain with a thickness of 180 m and a side of 2200 m (Figures 2 and 3). The simulated mine is placed just in the center of the model. Thus, the distance between the underground reservoir and the outer boundaries is of 1000 m. This distance allows to minimize the effects of the outer BCs on the results. The underground reservoir is represented by the nine underground chambers as described above (CH1 to CH9) that are hydraulically connected by galleries. The operation shaft is modeled using a rectangular prism adjacent to CH1 that connects the mine to the surface (Figures 2 and 3).

ures 2 and 3).

ters are typical of slate mines [22,37].

Thus, the distance between the underground reservoir and the outer boundaries is of 1000 m. This distance allows to minimize the effects of the outer BCs on the results. The underground reservoir is represented by the nine underground chambers as described above (CH1 to CH9) that are hydraulically connected by galleries. The operation shaft is modeled using a rectangular prism adjacent to CH1 that connects the mine to the surface (Fig-

The chambers and the operation shaft are modelled as independent domains that behave as linear reservoirs. Hydrogeologically, the groundwater flow behavior in the Martelange site is governed by fractures. However, it is known that a fracture-based groundwater flow is particularly difficult to implement in classical groundwater numerical models. In these types of cases, the underground medium is modeled using an equivalent porous medium (EPM) approach [42]. EPM approaches have been proven to be suitable to estimate the global groundwater behavior by numerous authors, including [48,49], particularly in sites with the presence of a high density of fractures as in the study site.

Concerning the hydraulic parameters, the *K* was modified in four different scenarios to observe its influence on the system behavior (note that the given EPM approach is considered, and *K* refers to the equivalent hydraulic conductivity). The values adopted for *K* in the different scenarios are specified in Section 2.2.3. The other hydraulic parameters

**Figure 3.** A detailed view of the modelled mine is shown on the right while a general view of the numerical model is shown on the left side of the figure. It is possible to observe how the bottom of the chambers increases progressively from west to east. **Figure 3.** A detailed view of the modelled mine is shown on the right while a general view of the numerical model is shown on the left side of the figure. It is possible to observe how the bottom of the chambers increases progressively from west to east.

The spatial and temporal discretizations were as follows: The domain was divided vertically in 29 layers and horizontally discretized using 3D prismatic elements. The maximum horizontal size of the elements was 150 m at the outer boundaries, and the mesh was refined around the mine where the horizontal size was about 5 m (Figure 3). The model was composed of 64,844 nodes and 38,680 elements. The total simulation length was one year and was divided into constant time steps of 15 minutes. Larger time steps would induce convergence problems and errors. Three different types of BCs were implemented. First, Dirichlet BCs were used to The chambers and the operation shaft are modelled as independent domains that behave as linear reservoirs. Hydrogeologically, the groundwater flow behavior in the Martelange site is governed by fractures. However, it is known that a fracture-based groundwater flow is particularly difficult to implement in classical groundwater numerical models. In these types of cases, the underground medium is modeled using an equivalent porous medium (EPM) approach [42]. EPM approaches have been proven to be suitable to estimate the global groundwater behavior by numerous authors, including [48,49], particularly in sites with the presence of a high density of fractures as in the study site.

prescribe the piezometric head on the W and E outer boundaries of the model. Two different hypotheses (scenarios TOP and MIDDLE) were considered to assess the influence of the elevation of the piezometric head on the system performance. In the TOP scenarios, the piezometric head was fixed at an elevation (i.e., with respect to the bottom of the model) of 121 and 120 m on the upgradient (W) and downgradient (E) sides, respectively. Consequently, the mine was totally full of water in the natural conditions. In the MIDDLE scenarios, groundwater head was prescribed at 76 and 75 m on the upgradient Concerning the hydraulic parameters, the *K* was modified in four different scenarios to observe its influence on the system behavior (note that the given EPM approach is considered, and *K* refers to the equivalent hydraulic conductivity). The values adopted for *K* in the different scenarios are specified in Section 2.2.3. The other hydraulic parameters were the same in all scenarios: the saturated and residual water contents were 0.05 and 0.01, respectively, while the specific storage coefficient was 10−<sup>4</sup> m−<sup>1</sup> . The chosen parameters are typical of slate mines [22,37].

(W) and downgradient (E) sides, respectively. In this case, only half of the volume of the mine was full of water under natural conditions. As a result of the prescribed heads, the The spatial and temporal discretizations were as follows: The domain was divided vertically in 29 layers and horizontally discretized using 3D prismatic elements. The maximum horizontal size of the elements was 150 m at the outer boundaries, and the mesh was refined around the mine where the horizontal size was about 5 m (Figure 3). The model was composed of 64,844 nodes and 38,680 elements. The total simulation length was one year and was divided into constant time steps of 15 minutes. Larger time steps would induce convergence problems and errors.

> Three different types of BCs were implemented. First, Dirichlet BCs were used to prescribe the piezometric head on the W and E outer boundaries of the model. Two different hypotheses (scenarios TOP and MIDDLE) were considered to assess the influence of the elevation of the piezometric head on the system performance. In the TOP scenarios, the piezometric head was fixed at an elevation (i.e., with respect to the bottom of the model) of 121 and 120 m on the upgradient (W) and downgradient (E) sides, respectively.

> Consequently, the mine was totally full of water in the natural conditions. In the MIDDLE scenarios, groundwater head was prescribed at 76 and 75 m on the upgradient (W) and downgradient (E) sides, respectively. In this case, only half of the volume of the mine was full of water under natural conditions. As a result of the prescribed heads, the hydraulic gradient was 4.6 <sup>×</sup> <sup>10</sup>−<sup>4</sup> for both hypotheses, and the groundwater flowed from W to E (Figure 2). The BCs implemented at the outer boundaries did not change through the simulations. Secondly, no-flow BCs were assigned to the top and the bottom of the modeled domain and to the N and S boundaries.

> Finally, pumping and discharge operations were simulated by use of a Neuman BC prescribing discharged or pumped water from the underground reservoir through the operation shaft. The value of the pumping and discharge rates were 5.94 m3/s, which is

the required flow rate to fill or empty 10% of the underground reservoir in 2 hours. The frequency of pumping and discharge (operation scenario) was generated randomly since it is difficult to predict how the electrical generation and demand will evolve during a year. Every two hours, a choice was made between three options (discharge, pumping, or no-operation), and thus, the minimum duration of pumping or discharge operations was 2 hours. No limitation was adopted concerning the duration of pumping, discharge, or no-operation phases. the required flow rate to fill or empty 10% of the underground reservoir in 2 hours. The frequency of pumping and discharge (operation scenario) was generated randomly since it is difficult to predict how the electrical generation and demand will evolve during a year. Every two hours, a choice was made between three options (discharge, pumping, or no-operation), and thus, the minimum duration of pumping or discharge operations was 2 hours. No limitation was adopted concerning the duration of pumping, discharge, or no-operation phases. Figure 4 displays the operation scenarios randomly computed for the TOP and MID-

hydraulic gradient was 4.6 × 10−4 for both hypotheses, and the groundwater flowed from W to E (Figure 2). The BCs implemented at the outer boundaries did not change through the simulations. Secondly, no-flow BCs were assigned to the top and the bottom of the

Finally, pumping and discharge operations were simulated by use of a Neuman BC prescribing discharged or pumped water from the underground reservoir through the operation shaft. The value of the pumping and discharge rates were 5.94 m3/s, which is

*Appl. Sci.* **2021**, *11*, x FOR PEER REVIEW 8 of 16

modeled domain and to the N and S boundaries.

Figure 4 displays the operation scenarios randomly computed for the TOP and MID-DLE hypotheses during the 10 first days and assuming no-groundwater exchanges. Positive values indicate that water was discharged while negative values indicate that the water was pumped. The same pumping–discharge function was used for all scenarios within the same hypothesis concerning the position of the piezometric head (i.e., TOP or MID-DLE). Later, during the simulation process, the virtual connections (the internal boundary conditions explained in Section 2.2.1) constrained the pumping and discharge when the underground reservoir was filled or emptied faster than expected. This occurred because water exchanges cannot be anticipated and, thus, cannot be taken into account when defining the operation scenarios. DLE hypotheses during the 10 first days and assuming no-groundwater exchanges. Positive values indicate that water was discharged while negative values indicate that the water was pumped. The same pumping–discharge function was used for all scenarios within the same hypothesis concerning the position of the piezometric head (i.e., TOP or MID-DLE). Later, during the simulation process, the virtual connections (the internal boundary conditions explained in Section 2.2.1) constrained the pumping and discharge when the underground reservoir was filled or emptied faster than expected. This occurred because water exchanges cannot be anticipated and, thus, cannot be taken into account when defining the operation scenarios.

**Figure 4.** Operation scenarios during the first 10 days for the TOP (blue continuous line) and MID-DLE (red dotted line) scenarios. **Figure 4.** Operation scenarios during the first 10 days for the TOP (blue continuous line) and MIDDLE (red dotted line) scenarios.

Figure 4 shows that the pumping–discharge frequencies were very similar for the TOP and MIDDLE scenarios. Differences were observed only during the first and third days. In the first day, we observed that the discharged water was higher for the MIDDLE scenario (red dotted line), and this occurred due to the BCs and the initial conditions. Initially, the underground reservoir was full in the scenario TOP, and, therefore, pumping was needed before any discharge. Figure 4 shows that the pumping–discharge frequencies were very similar for the TOP and MIDDLE scenarios. Differences were observed only during the first and third days. In the first day, we observed that the discharged water was higher for the MIDDLE scenario (red dotted line), and this occurred due to the BCs and the initial conditions. Initially, the underground reservoir was full in the scenario TOP, and, therefore, pumping was needed before any discharge.

There was a short pumping at the beginning, and then the underground reservoir was quickly filled by a short subsequent discharge; therefore, no more water could be There was a short pumping at the beginning, and then the underground reservoir was quickly filled by a short subsequent discharge; therefore, no more water could be introduced in the underground reservoir until the next pumping. However, in the scenario MIDDLE, the initial head was located at the half depth of the mine; therefore, more water could be discharged without the need for previous pumping. In contrast, during the third day, more water could be pumped in the scenario TOP than in the scenario MIDDLE. In this case, the natural piezometric head in the scenario MIDDLE was lower than in the TOP one; therefore, the underground reservoir was empty faster in the MIDDLE than in the TOP scenario, and the pumping was stopped earlier.

For the initial conditions, we assumed that the piezometric head distribution matched with that in natural conditions. Therefore, no previous pumping was considered before the start of the activity of the plant. Certain authors considered that, initially, the underground reservoir must be empty [50] because previous works needed to adapt the mine to be used as an underground reservoir, and this adaptation requires dewatering it. However, it was proven in a previous investigation that the differences between considering a previous pumping or not are negligible and only occur during the early periods [33].

This justifies the adopted hypothesis for constructing the numerical model used here. In addition, after any long shutdown of the plant activity, the piezometric head would reach its natural position. As a result, the initial conditions when the activity of the plant would be resumed are the same as those considered in the numerical model. It is important to clarify that the initial saturated thickness varies depending on the simulated scenario (TOP or MIDDLE), as they have different BCs on the outer boundaries.
