## 2 Materials and Methods
### Analytical Model
\[\frac{dT_{1}}{dt}=\frac{2U_{h}}{m_{s}C_{ps}}(T_{a}-T_{1})-\frac{2U_{s}}{m_{s} C_{ps}}(T_{1}-T_{2})-\frac{dT_{2}}{dt} \tag{5}\]  
\[\frac{dT_{2}}{dt}=\frac{2U_{s}}{m_{c}C_{pc}}(T_{1}-T_{2})-\frac{2U_{c}}{m_{c}C _{pc}}(T_{2}-T_{3}) \tag{6}\]  
\[\frac{dT_{3}}{dt}=\frac{2U_{c}}{m_{ss}C_{pss}}(T_{2}-T_{3})-\frac{2U_{ss}}{m_{ ss}C_{pss}}(T_{3}-T_{4}) \tag{7}\]  
where \(U_{h}\), \(U_{s}\), \(U_{c}\), and \(U_{ss}\) are the convection transmittance, sealing layer transmittance, concrete transmittance, and sandstone transmittance, respectively. These thermal coefficients are evaluated according to  
\[U_{h}=h2\pi r_{1}L \tag{8}\]  
\[U_{s}=\frac{2\pi K_{s}L}{Ln\bigg{(}\frac{r_{2}}{r_{1}}\bigg{)}} \tag{9}\]  
Figure 2:0  
Figure 2:1  
In these expressions, \(h\) is the film coefficient of heat transfer, in Wm\({}^{-2}\) K\({}^{-1}\); \(L\) is the tunnel length, in m; \(K_{s}\), \(K_{c}\) and \(K_{ss}\) are the thermal conductivity of sealing layer, concrete and sandstone, in Wm\({}^{-1}\) K\({}^{-1}\); \(C_{ps}\), \(C_{pc}\), and \(C_{pss}\) are the specific heat at constant pressure for the sealing layer, reinforced concrete and sandstone, in J kg\({}^{-1}\) K\({}^{-1}\); and \(r_{1}\), \(r_{2}\), \(r_{3}\), and \(r_{4}\) are the (equivalent) radius of sealing layer, concrete, sandstone and external walls, in m (Figure 2b). Note that due to the geometrical characteristics of the tunnel, the heat transfer formulation in cylindrical coordinates has been used for Equations (8)-(11). For the evaluation of the convection film coefficient, the correlation proposed by [PERSON] et al. [32] and [PERSON] et al. [33] for the measurement of heat transfer coefficients in high-pressure vessels during gas charging was considered. In particular, a mixed (natural and forced) convection is modelled as a combination of both unsteady Reynolds and Rayleigh numbers. For the present investigation, a slight correction has been introducedfor the exponent of the Rayleigh number, resulting the Nusselt number in the following expression:  
Figure 2:2  
where \(Re\) is defined as the Reynolds number of the incoming mass flow rate (\(Re=4\mu/2\mu\pi r_{1}\)) and \(Ra\) is computed from the instantaneous thermal properties of the air (\(Ra=g\beta(T_{a}-T_{1})C_{p}\rho^{2}L^{3}/(\mu K)\)), being \(g\) the gravity acceleration, \(\beta\) is the volumetric thermal expansion coefficient, \(Cp\) the specific heat coefficient at constant pressure, \(\rho\) the density, \(L\) the tunnel length, \(\mu\) the dynamic viscosity and \(K\) the thermal conductivity. The film coefficient is thus computed as  
Figure 2:3  
Finally, a backward Euler explicit discretization has been employed over Equations (3)-(7) to resolve the coupled system of equations. A typical time-step of 0.01 s was applied to provide an accurate temporal description of the different temperatures. The model was resolved iteratively until a maximum prescribe pressure was reached inside the reservoir. In addition, according to [PERSON] et al. (2018) and [PERSON] et al. (2015), the air has been considered as a real gas through a compressibility factor (Z) that it is computed using the Berthelot gas state equation. Finally, considering a density value uniformly distributed in the volume, the pressure value of the CAES may be obtained at any instant by applying Equations (14) and (15). Figure 2:4  
Figure 2:5  
where \(T_{c}\) and \(P_{c}\) the air temperature and the air pressure at critical conditions, assumed to be 132.65 K and 3.76 MPa. The contact surfaces between air-sealing layer (\(A_{1}=2\pi r_{1}L\)), sealing layer-concrete (\(A_{2}\)-\(2\pi r_{2}L\)), concrete-sandstone (\(A_{3}\)-\(2\pi r_{3}L\)), and sandstone-exterior (\(A_{4}\)-\(2\pi r_{4}L\)), are 359, 365, 545, and 1610 m\({}^{2}\), respectively.