## 2 Materials and Methods
### Methods
#### 2.3.2 Conceptualization  
A total of 20 layers were discretized in the 3D model. From the top, ten layers represent the shallow confined Cenozoic aquifer system: the first 5 layers correspond to the top part of the upper aquifer, layers 6 and 7 are used for a clayey part, and layers 8-10 for the lower aquifer sandy part of these Cenozoic sediments. Three layers are used to represent the aquitard consisting of the clay layer at the bottom of the Cenozoic sediments, the Cretaceous marls, and the Cretaceous chalks of low hydraulic conductivity. Then, seven layers represent the fractured Paleozoic bedrock until a depth of 150 m: five layers in the upper fractured zone, and two layers in the less fractured and less permeable zone. The bottom of the model is arbitrarily set at an altitude of \(-\)150 m (i.e., a depth of about 170 m) where relatively low permeability is observed (by pumping tests) and because the screens of the deep ATES system (Building 3) are located at an altitude between \(-\)67 m and \(-\)138 m (Figure 3). The groundwater flow is simulated in 3D including vertical and horizontal flow components (cfr Section 3.1). Feflow(r) allows to simulate the effect of temperature changes on the groundwater density and viscosity. However, constant hydraulic conductivity values have been adopted as the variation can be neglected for limited changes in temperature.  
The boundaries and boundary conditions (BCs) were chosen according to regional hydrogeological data [33]. The top of the model corresponds to a no-flow BC, which is induced by the low permeability of the overlying confining layers and the largely urbanized land surface. As mentioned previously, the bottom boundary is set at a depth of 150 m with a no-flow BC as at this depth, the bedrock is less permeable [34] and the vertical component of the flow can be neglected. The lateral boundaries were chosen at a 'large enough' distance compared with the radius of influence as interpreted from local pumping tests. On the SE and NW lateral boundaries (Figure 4), potentiometric heads are prescribed along head contour lines in the shallow aquifer and based on interpolated heads from the potentiometric map in the bedrock aquifer [33]. No flux BCs were chosen on the NE and SW lateral boundaries (Figure 4) as they are globally parallel to the general groundwater flow directions, as well in the shallow aquifer and also in the deep aquifer. All the ATES pumping and injection wells and other water production wells are represented by prescribed fluxes (i.e., depth-averaged distributed onto different layers if different layers are screened in the borehole) on the corresponding nodes of the finite element mesh.  
Figure 3: 3D view of the finite element mesh and a detailed NW–SE cross-section passing through the location of Building 3 (deep ATES system).  
In the development zone (in yellow in Figure 4), the size of the finite element can be very small, from a few cm near the wells to a few tens of meters near the boundaries of this area of interest. In total, the 3D model consists of 1,248,420 finite elements with 20 layers of 62,421 elements. Specifically, in the development zone (area of interest), the finite element mesh refinement leads to a total number of finite elements of 434,520.  
For the heat transport BCs, a temperature of 12.4 \({}^{\circ}\)C was chosen along the lateral boundaries. This temperature is the natural background temperature. A zero heat flux is assumed on the top as on the bottom of the model. This means that the air temperature influence and the mean deep geothermal energy flux are neglected [3; 35].  
As a first step, a steady-state groundwater flow was calculated from the initial potentiometric conditions. Then, the calculated steady-state potentiometric map was used for starting the transient groundwater flow. For heat transport, the natural background temperature was chosen as the initial temperature in the whole model.