Joule-Thomson cooling of CO2 injected into aquifer under heat exchange with adjacent formations by Newtons law- 1D exact solution

This paper discusses axi-symmetric flow during CO2 injection into a non-adiabatic reservoir accounting for Joule-Thomson cooling and steady-state heat exchange between the reservoir and the adjacent layers by Newtons law. An exact solution for this 1D problem is derived and a new method for model validation by comparison with quasi 2D analytical heat-conductivity solution is developed. The temperature profile obtained by the analytical solution shows a temperature decrease to a minimum value, followed by a sharp increase to initial reservoir temperature. The analytical model exhibits stabilization of the temperature profile and cooled zone. When accounting for heat gain from adjacent layers, the minimum temperature is higher compared to the case with no heat exchange. The analytical model can be used to evaluate the risks of hydrate formation and/or rock integrity during CO2 storage in depleted reservoirs.


𝑎 𝑎
Thermal diffusivity of adjacent formations m 2 /s

INTRODUCTION
The growing concerns over negative impacts of greenhouse gases, in particular Carbon dioxide (CO2), on climate change have increased the interest in developing different methods to reduce their emissions into the atmosphere.
A practical and complementary solution is to capture CO2 from highly concentrated sources, e.g., industrial plants or directly from air and store it in underground geological formations such as depleted oil and gas reservoirs, deep saline aquifers, and coal beds.
Depleted gas fields exhibit several advantages over other geological formations.The extracted gas from these fields can in principle be replaced by CO2 and as such, huge volumes of CO2 can be stored in these reservoirs.For example, it has been estimated that in the Dutch sector of the North Sea, theoretically, more than 1.5 Gton of CO2 can be stored in the depleted gas fields (Van der Velde et al., 2008).Due to large recovery factors (extracted fraction of the initial volumes) and compressibility of gas more space is available in the gas fields compared to oil fields, in which the extracted oil is usually replaced by injection of another (usually incompressible) fluids (Hamza et al., 2021).
The production of oil fields often involves injection of external fluids (for pressure maintenance or improvedrecovery purposes) and/or drilling of a much larger number of wells, limiting their storage volumes and increasing the risk of leakage pathways (Van der Velde et al., 2008).Compared to aquifers, because of their long production history, the uncertainty in the geological settings (permeability, heterogeneity, faults, fractures, etc.) of the hydrocarbon fields is relatively low.Moreover, these fields have proven seal and containment integrity and have part of the infrastructure required to handle the gas.The sequestration of CO2 can also be combined with enhancing gas recovery, which could partially compensate for the cost associated with CO2 sequestration (Battashi et al., 2022;Hamza et al., 2021).
However, there are challenges with CO2 storage in depleted gas reservoirs, some of which are related to the thermodynamic properties of CO2.From point of capture (or production) to inside the reservoir pores, CO2 experiences different pressure and temperature conditions and may exist in gas, liquid or, under certain conditions, solid forms.Of particular concern for CO2 storage in low-pressure reservoirs is the so-called Joule-Thomson (J-T) cooling effect, which is caused by expansion of CO2 from high surface pressures to low reservoir pressures.This is schematically shown by path A→D in Fig. 1.Point A represents the pressure and the temperature of the transported CO2 at the storage site, which is assumed to be at 100 bar and 30 o C (i.e., CO2 is injected in dense liquid phase).With the assumptions that heat conduction in the wellbore is negligible and that the viscous and gravity forces balance each other, the bottom-hole pressure of CO2 is assumed to be close to that of point A (although in practice, the downhole temperature of CO2 will be slightly higher due to heat conduction from the well, and CO2 will most likely be in two phase regime).Isenthalpic expansion of CO2 from 100 bar to a reservoir pressure of 20 bar (point D) reduces the CO2 temperature to -10 o C. In general, for reservoirs with pressures lower than 35 bar J-T cooling effect reduces the CO2 temperature to sub-zero values, as illustrated in Fig. 1.This can potentially lead to the formation of hydrates and freezing of pore water, as implied from the pressure-temperature diagram of the water/CO2 system in Fig. 2a.The near-wellbore appearance of hydrates or ice can in turn result in the impairment of well injectivity and/or loss of well-related containment (Vilarrasa & Rutqvist, 2017).Moreover, the thermal stress generated by the cooling effect can lead to fracturing of the rock, the extent of which needs further investigation.
An alternative solution will be to heat CO2 before injection (path A→B→C in Fig. 1); however, this method is energy (and CO2) intensive and adds further to the cost of CO2 storage.Injection of CO2 in gas form is also not desirable because storing similar mass rates of CO2 require the drilling of multiple wells and requires a larger pore space.The long-term and safe storage of CO2 requires a detailed understanding of the physical, chemical, geo-mechanical and thermal effects caused by the injection of CO2 in the reservoirs.This includes mobilisation of fines by CO2water menisci (Chequer et al., 2020) , salinity alteration (Othman et al., 2019;Parvan et al., 2021;Yu et al., 2018) hydrate formation (Machado et al., 2023) and salt precipitation (Miri et al., 2015;Moghadasi et al., 2004).In addition, to identify the type of risk associated with cold CO2 injection and design a mitigation plan, it is crucial to quantify the range of expected temperatures in the depleted gas reservoirs.Oldenburg (2007) found that for a constant injection rate, lower permeability and higher porosity increase the effect of J-T cooling.Mathias et al. (2010) derived an exact solution for CO2 injection accounting for J-T cooling into an adiabatic reservoir (ignoring heat exchange with surrounding formations).The analytical model presented in his paper provides explicit formulae for the temperature and pressure profiles, and position of the cold front and quantifies the impact of several parameters on the cooling effect.This model, like all analytical models, allows for fast calculations for multivariant sensitivity studies and can benchmark the numerical methods for more complex modelling (Kacimov & Obnosov, 2023a, 2023b;Katzourakis & Chrysikopoulos, 2019;Moreno et al., 2021).These advantages explain numerous studies on the analytical modelling of problems related to CO2 storage (Ahmadi & Chen, 2019;Celia et al., 2011;Moreno et al., 2021;Ziv Moreno & Avinoam Rabinovich, 2021;Z. Moreno & A. Rabinovich, 2021;Nguyen et al., 2022;Norouzi et al., 2022).
Heat exchange between the reservoir and the adjacent layers can significantly affect the non-isothermal flow in porous media (Bedrikovetsky, 1993;Lawal, 2020).Nevertheless, an exact solution for CO2 injection into lowpressure reservoirs or aquifer accounting for this effect is not available.This paper drives an exact solution accounting for heat exchange between the reservoir and adjacent layers.
The explicit formula allows studying the timely evolution of the temperature profile, including the temperature penetration depth.The solution shows that the temperature profile stabilises with time and the penetration front stops.The solution also allows for calculating the well pressure at the injection well or well injectivity index and the dynamics of the heat penetration depth The explicit formulae for temperature and pressure allow calculating maximum injection rate that avoid formation of hydrates, i.e. the path well-reservoir in phase diagram does not enter the hydrate domain (Fig. 2).The validity of the model with Newton's law for heat exchange is determined by comparison with the exact solution.
The structure of the paper is as follows.Section 2 describes the major model assumptions and presents the mathematical model and the ensuing equations.Section 3 presents the explicit solution of the model, derived in Appendix A. Section 4 provides the calculation for the area of validity of the steady-state heat exchange model.
Appendix B presents the exact solution for heat flux in adjacent layers and derives the criterion used to establish the domain of validity for the steady-state heat exchange term used in this paper.Section 5 presents the results of the analytical model.Section 6 discusses the sensitivity study of the model, introducing the impact of dimensionless parameters on the temperature.

MATHEMATICAL MODEL
This section presents the mathematical model for temperature profile evolution during CO2 injection: the model assumptions (section 2.1) and derivation of the heat-transfer equation (section 2.2).

Model Assumptions
The main assumptions of the model are:  The model assumes steady-state heat exchange between the reservoir and the surrounding formations (Newton's law) 2   −1 ( −   ) (Jang et al., 2022).So, the reservoir is bounded by overburden and underburden shales with thickness  [L].Here, r [L] is the radial distance from the injection well, (, ) [K] is the fluid temperature in the reservoir,   [MLT -3 K -1 ] is the heat-transfer coefficient (or overall thermal conductivity) of the seal in contact with CO2, and TI [K] is the initial temperature of the reservoir.This assumption is valid when the reservoir is hydrodynamically separated above and below from the adjacent layers by the impermeable seals.The heat fluxes behind and ahead of the horizontal boundary between the adjacent formations and seals are assumed equal.Heat flux is a product of thermal conductivity and temperature gradient.The thermal conductivity of an impermeable seal, e.g., shale is determined by that of the solid matrix, which is significantly lower than the thermal conductivity of the water residing in the adjacent permeable formation and forming a conductive cluster.Therefore, the temperature gradient in shales is expected to be significantly lower than the temperature gradient in the adjacent formation.This supports the assumption that the temperature at the upper and lower bounds of the adjacent shales layers remains equal to the initial temperature TI.
Newton's law for the heat exchange between the reservoir and the adjacent layers is a commonly-used assumption for 1D non-isothermal flows in porous media (Atkinson & Ramey, 1977;Batycky & Brenner, 1997;Fedorov & Sharafutdinov, 1989;Gordeev et al., 1987;Jang et al., 2022;LaForce et al., 2014;Muradov & Davies, 2012;Muradov & Davies, 2009;Payne & Straughan, 1998;Pires et al., 2006;Xu et al., 2013;Yortsos & Gavalas, 1982;Zazovskii, 1983;Zolotukhin, 1979).However, the physical-geological conditions for Newton's law validity are not present in the literature.Derivations in Appendix B present the dimensionless criteria determining the validity of Newton's law.Furthermore, it is assumed that the heat exchange between the injected CO2, water, and the rock occurs instantly.The heat transfer due to conduction in the flow direction is ignored compared to advection.The schematic of the radial problem for CO2 injection into a depleted reservoir is given in Fig. 3.

1-D Flow Equations
The energy-balance equation consists of the accumulation term of heat in the rock, CO2, and water, advective heat transport, the J-T effect, and heat exchange of the reservoir with the adjacent layers (Lake et al. 2014;Wang 2022; Hashish & Zeidouni 2022): where is the fluid pressure.The unknown in Eq. ( 1) is the reservoir temperature T(r,t).
The J-T coefficient of CO2,   , depends on the pressure and temperature as shown in Fig. 4. For the depleted reservoirs with pressures of less than 5 MPa,   increases with a decreasing temperature.For a fixed temperature,   can be assumed independent of pressure for these reservoirs, even though for larger pressures   decreases with increasing pressure.
The overall heat transfer coefficient  [MT -3 K -1 ] is defined as (Lawal, 2020):  2021) ) It has been experimentally observed that when a fluid with a different temperature than the temperature of the porous medium is injected the heat loss (or gain) per unit length is almost constant.Consequently, the overall heat transfer coefficient, , is assumed constant in this study.In other words, the heat conductivity in the z-direction is assumed to be very large, which results in a constant temperature across the vertical thickness of the reservoir (LaForce et al., 2014).However,  may vary with injection rate, the surface area in contact with surroundings, time, thermal conductivity of the surrounding porous medium, and/or the heat capacity of the reservoir, and the type of process (cold or hot fluid injection) (LaForce et al., 2014).Cold CO2 injection into a reservoir filled with a hot fluid is analogous to cold water injection in geothermal reservoirs.The experimental data show that for cold water injection lower injection rates are more advantageous for heat extraction from the surrounding rocks.
The volumetric injection rate,   , is defined as: where   [MT -1 ] is the mass injection rate.The position of the CO2 front   [L] is calculated assuming gas incompressibility: (5) The pressure gradient in the reservoir is calculated using Darcy's law: where   [ML -1 T -1 ] is viscosity of CO2, and    is the end-point relative permeability of CO2.Finally, the substitution of Eq. ( 6) into Eq.( 1) yields the first order partial differential equation: The boundary condition (BC) correspond to the injection of CO2 with a constant pressure and temperature: where   [L] is the well radius.
The initial condition (IC) corresponds to the equality of the initial temperatures of the main reservoir rock and the adjacent layers, i.e., The following dimensionless variables are introduced to nondimensionalize Eq. ( 7).
Further, the dimensionless constants V* and A* are defined as: By substituting Eqs.(10-12), the dimensional equation Eq. ( 7) transforms into the dimensionless form: Additionally, the BC Eq. ( 8) and IC Eq. ( 9) are also expressed in dimensionless form as: It is assumed that the heat transfer within the rock occurs instantaneously, i.e., the injected cold CO2 is instantly heated up by the residing water and rock.This implies that the temperature front lags significantly behind the CO2 front inside the reservoir.

EXACT ANALYTICAL SOLUTION
The exact solution of Eq. ( 13) with initial and boundary conditions Eqs. (14-15) is obtained by the method of characteristics (Boyce & DiPrima, 2013;Polyanin & Zaitsev, 2003)  Here, the dimensional constants V and B are: The analytical model (A16), ( 16) can be extended for the case where the injected temperature changes with time -TJ(tD).In this case, TJ in Eq. (A16) must be substituted by TJ(tD-(rD 2 -1)/2V*).
Let us define the depth of temperature wave penetration as a coordinate of the "centre" of the wave.For variable r, the probability distribution function can be defined as: The arithmetic mean of r is

AREA OF VALIDITY OF THE MODEL
This section determines the area of validity of the analytical model, including comparison between the solutions for adiabatic and non-adiabatic reservoirs (section 4.1), formulation of the method of 1D model validation by comparison with 2D solution (section 4.2), and calculations of validity intervals for the model parameters and independent variables (section 4.3).

Definition of the Basic Case for Sensitivity and Model Validation
To validate the developed analytical model, the solution given by Eq. ( 19) was compared with the results of Mathias et al. (2010) in the limit of no heat transfer between the reservoir and the adjacent layers, i.e.,  → 0. The parameters provided in Table 1 are used as a basic case for all calculations in this paper.Table 2 presents the intervals of variation of variables.Fig. 5 shows the close agreement between the numerical implementations of two models.
As expected, the temperature front falls behind the CO2 front because of the heat provided by the rock (the volumetric heat capacity of the rock is much larger than those of both water and CO2).The temperature decreases from injection temperature to a minimum value and then sharply increases to the initial temperature of the reservoir.
There are two distinctive features for the solution of the case, where κ →0: (1) temperature decreases to a minimum value just upstream of the sharp front at each timestep, and (2) the temperature front keeps moving inside the reservoir for as long as CO2 injection continues.The sharp increase to the initial temperature (shock front) is due to absence of conduction in the analytical model.In accordance with the properties of the radial flow (i.e., large pressure drops near the injection point), the cold front penetrates quickly inside the reservoir.If the injected CO2 is undersaturated, the residing water in the pores evaporates; and therefore, the water saturation reduces to near-zero close to the injection well.However, the speed of the temperature front is expected to be larger than that of the drying front such that a region with a non-zero water saturation appears between the dry-out and cold fronts in a relatively short time after injection of CO2.This is schematically shown in Fig. 6.If the pressure and temperature fall within in the hydrate region on the CO2/water phase diagram in Fig. 2, CO2 hydrates might form potentially leading to blockage of CO2 flow.The primary assumption of the model given by Eq. ( 1) is steady state heat exchange between surrounding formations and impermeable shale that confines the reservoir, which is the Newton's law.This section establishes the domain of validity for the heat exchange term used in the energy balance (Eq.( 1)), aiming to determine the condition under which the discrepancy between the temperature at the interface, (  = 0,   ), and the initial reservoir temperature,  =   →  = 1, is considered insignificant.In our case, we constrain the discrepancy between the two solutions for temperature to the difference of less than or equal to 10%.
We utilize the solutions obtained in Appendix B for the temperature profile at the interface between adjacent impermeable shale, as derived in Eq. (B19).The condition in dimensionless variables is as follows: In this analysis we introduce two dimensionless constants, A* and V* given in Eq. ( 11) and Eq. ( 12).Also, we present the constant ω determined in Appendix B as Eq.(B9): Newton's law is valid when the temperature on the outer bounds of the adjacent seals is equal to the initial reservoir temperature.This condition depends on conductivity and capacity of the adjacent formation layers i.e., depends on the heat conductivity (  ) and volumetric heat capacity (  ).Indeed, equality of the temperature on the outer shale bounds to initial reservoir temperature corresponds to ω tending to zero.
The intervals for parameters A*, V* and ω are determined from typical fluid and rock properties.The details below highlight the values adopted in this work.The heat capacities of gas, water and rock varying from 709 J (kg K) -1 to 1476 J (kg K) -1 , 3965 J (kg K) -1 to 4335 J (kg K) -1 and from 776 J (kg K) -1 to 1215 J (kg K) -1 , respectively.The heat conductivities of shale and adjacent formations ranges from 1 W (m K) -1 to 3.7 W (m K) -1 and from 3.7 W (m K) -1 to 5 W (m K) -1 , respectively.The densities of gas, water and rock vary from 0.5 kg m -3 to 1236 kg m -3 , from 527 kg m -3 to 1000 kg m -3 and from 2270 kg m -3 to 3200 kg m -3 , respectively.In addition, the thickness of the reservoir and the shale vary from 1 m to 100 m and from 1 m to 3.33 m, respectively.Well radius varies from 0.05 m to 0.13 m.
Given that our solution in Eq. ( B19) is notably influenced by system parameters-A*, V* and ω, a domain made of the three dimensionless constants, (A*, V*, ω) in Eq. ( B19) is established.The eight vertices of the tetrahedron represented by dots in Fig. 7 correspond to the extreme values of A*, V* and ω.The inner dot numbered 9 corresponds to the unperturbed base case given in Table 1Error!Reference source not found.. Table 2 summarizes the critical values of the 9 coordinate points illustrated in Fig. 7. Our aim is to determine a region of validity where the condition in Eq. ( 21) is upheld.To visualize this, we plot   1 − (  = 0,   ) < 0.1.In case a different criterion is needed, the area of validity for Newton's law changes accordingly.To provide an analysis on the validity in each of the points presented in Table 2, we introduce the definition for radius of injection, following Eq.( 5):

RESULTS OF THE ANALYTICAL MODELLING
This section presents the results of the derived analytical model and discusses their relevance for CO2 storage in the low-pressure or depleted fields.The impact of heat exchange between the main formation and the surrounding rock on the propagation of the temperature front is shown in Fig. 11.A heat transfer coefficient of κ = 2.5 W/m 2 /K was used in the calculations (Cermak & Rybach, 1982;Labus & Labus, 2018;Robertson, 1988;Schoen, 2015).For the purpose of comparison, Fig. 11 further shows the case of no heat exchange ( → 0) with the dashed curves.Initially and for relatively short times, the impact of the heat gain from the surrounding layers is not significant.
Fig. 11.The impact of heat exchange on temperature profile (solid lines) in comparison to adiabatic system (dashed lines), CO2 front position and Pressure drop for variying times.
However, as CO2 injection continues, the heat supplied by these layers heats the injected CO2 and therefore the temperature increases slightly before jumping to the initial temperature of the reservoir.Unlike the case with no heat exchange, the temperature upstream of the shock is not the lowest temperature.Indeed, a key difference is that the position and the value of the minimum temperature remains fixed for all times in the new model.The area between the solid and the dashed curves in Fig. 11 represents the extent of the impact of the heat exchange.The area increases with time and as a result the length of the shock front becomes smaller every timestep, until it disappears completely.This corresponds to the time at which the system reaches steady state, and the temperature front no longer moves.The solution of the ordinary differential equation (ODE) is obtained by the separation of variables: Here, C and D are defined as follows: The time intervals in Fig. 11 correspond to dimensionless timey, tD = 0.01, 0.09, 0.3, 1.8, and 5.4, respectively.Eqs.
(24) and ( 25) indicate that the position of the steady state temperature depends largely on the injection rate, the volumetric heat capacity of CO2 and rock properties.The pressure profile under the assumptions of gas incompressibility and constant viscosity is steady state; it is obtained by integration of Eq. ( 6) with BC Eq. ( 8) for pressure.The temperature profile in Fig. 11 shows temperature decreases from the well to the temperature front, and then monotonically increases until initial temperature.Minimum temperature is achieved exactly behind the temperature front.It is crucial to determine a minimum injection rate at which the temperature behind the temperature front and pressure along the temperature front are not conducive for hydrate formation (as illustrated by the phase diagram in Fig. 2a).Additionally, Fig. 2b.shows the pressure-temperature path in the phase diagram, where the reservoir point corresponds to minimum temperature.The length of the segment, injection well-reservoir, as it follows from Eq. ( 6), is proportional to injection rate qJ.This allows determining the minimum rate of hydrate formation.Fig. 2c shows two solutions that corresponds to mass injection rates of   = 3 kg/s and   = 4 kg/s.
According to Fig. 2c, hydrates formation occurs when the plot crosses the Hydrates + Gas CO2 phase boundary, which is evident only for the higher injection rate.

SENSITIVITY STUDY OF THE MODEL
This section investigates the sensitivity of the analytical model to the model parameters: heat transfer coefficient (subsection 6.1), permeability (subsection 6.2), and reservoir thickness as well as injection rate (subsection 6.3).
Additionally, the effect of the dimensionless parameters A* and V* on the temperature profile is investigated in subsection 6.4.

Effect of Heat Transfer Coefficient
The impact of the heat transfer coefficient, , is shown in Fig. 12.With the increase of , a higher amount CO2 is heated and therefore the cold front is retarded, compared to the case where  is smaller.Moreover, the minimum temperature is slightly lower for the case with the smaller , as illustrated in Fig. 12a.Fig. 12b shows that in the adiabatic reservoir, the penetration depth, which was introduced by Eq. ( 20), does not stabilize (blue curve).Stabilization of the heat penetration is achieved with the presence of heat exchange in the system.For the given temporal range, the curve corresponding to the lower value for heat conductivity does not stabilize, while the curve corresponding to the higher value achieves stabilization.The three values of heat conductivity presented in Fig. 12 correspond to the following dimensionless parameters-(A*,V*→∞), (A*=7.2×10 2 , V*=8.6×10 4 ), (A*=1.4×10 2 , V*=1.7×10 4 ), respectively.
6.2 Effect of Permeability Fig. 13 shows the sensitivity of the model to varying rock permeabilities.As stated earlier, the magnitude of temperature reduction due to the J-T effect depends on the pressure gradient and the J-T coefficient.For fixed injection pressure and flowrate at the boundary, the pressure gradient increases with decreasing reservoir permeability, causing an increased temperature drop.However, this conclusion is valid only when equal flowrates are considered for different permeabilities.For a heterogeneous reservoir, the injected CO2 will be distributed in proportion to the permeability of each layer.If the permeability contrast is large, then the large fraction of CO2 will flow through the high permeability layer causing a larger pressure drop compared to the low-permeability layer.
Therefore, for heterogeneous reservoirs with large permeability contrast, the "isenthalpic" expansion of CO2 will result in lower temperatures in the high permeability layers.
It is worth mentioning that when CO2 becomes liquid or the pressure exceeds the critical pressure of CO2 the magnitude of   also reduces, especially in the colder regions near wellbore (see Fig. 4) resulting in lower cooling effects.Furthermore, the presence of impurities in the injected CO2 or contamination of CO2 with the gas already present in the reservoir will also affect the temperature drop and the extent of the cooling zone (Ziabakhsh-Ganji & Kooi, 2014).Gases like SO2 expand the cooling zone while gases like N2 and CH4 tend to contract this zone (Ziabakhsh-Ganji & Kooi, 2014).The Joule-Thomson cooling effect is more pronounced at low-permeability and thin layers, assuming a constant injection rate.With the increase of the reservoir permeability or thickness the impact becomes less and after a certain value the steady-state temperature profile appears to be independent of these parameters.
Fig. 13.The impact of permeability on the temperature profiles.

Effect of System Parameters on Minimum Reservoir Temperature
In the sensitivity analysis shown in Fig. 14, only the investigated model parameter was changed while the other parameters remained constant.For simplicity the parameters are normalized to the parameters of the base-case presented in Table 1.Moreover, the rnorm on the y-axis is relative to the position of the steady-state temperature of the base case.The normalized distance has a near linear relationship with the mass injection rate, i.e., the higher the injection rate the longer the cold front travels inside the reservoir.In addition, the minimum temperature significantly decreases as the injection rate increases.For example, with doubling of   , the minimum temperature drops from 10 o C to -10 o C, which has huge implications for a CO2 storage project.

Sensitivity to Change in the Dimensionless Constants
The dimensionless constants A* and V*, as defined in Eq. ( 11) and Eq. ( 12), are the parameters that account for the J-T effect, and the heat front velocity, respectively.An increase in parameter A* is equivalent to an increase in the J-T coefficient, injection rate and shale thickness, as well as a decrease in reservoir thickness, shale heat conductivity and permeability.Also, an increase in parameter V* corresponds to an increase in shale thickness and injection rate, as well as a decrease in shale conductivity.Fig. 15 show the sensitivity of the temperature profile to the dimensionless constants, A* and V*.From Fig. 15a, it is observed that, while maintaining all other system parameters constant and varying A*, the temperature drop in the system ranges from ΔT = 170°C to ΔT = 0.In Fig. 15b, varying the values of V*, while holding all other system parameters fixed results to temperature drops varies between ΔT = 95°C to ΔT = 0. Hence, the system is more sensitive to modifications in parameter A* when assessing the temperature profiles.

DISCUSSION
The analytical model presented in this paper (Eq.( 19)) assumes a constant injection rate; under given pressure drawdown between the injector and the reservoir.Eq. ( 17) defines the dimensional parameter for heat front velocity.
As temperature front moves significantly slower than the CO2-water front; V ≪ 1, for which water saturation in the cold region is only slightly above Swi.However, varying two-phase mobility during the displacement can affect pressure behaviour, so well injectivity prediction requires accounting for the overall saturation profile.Self-similar solution for displacement of water by CO2 can be matched with the analytical model (Eq.( 19)).
Another mechanism requiring two-phase formulation during CO2 injection, which is not accounted for in this study, is mobilisation and migration of natural reservoir fines, that are detached by gas-water menisci (Chequer et al., 2020;Nguyen et al., 2022).The governing system includes the equation for motion of gas-water interface, where fines detachment occurs (Shapiro, 2015).The detachment is determined by the interaction of capillary and DLVO forces exerting the attached fine particle (Yuan & Moghanloo, 2018;Yuan & Moghanloo, 2019).The saturation variation from 1 to Swi around CO2-water front requires generalisation of the overall unsteady-state model for two-phase flow with moving interface (Shapiro, 2016(Shapiro, , 2018)).Accounting for two-phase flow with fines migration in the model Eq. ( 19) captures the effects of injectivity decline during CO2 injection.
Other mechanisms contributing to formation damage and injectivity decline, triggered by J-T cooling effect, encompass salt precipitation, rock dissolution, mineral precipitation reaction, and hydrate formation (G Moghanloo et al., 2017;Ge et al., 2020;Nguyen et al., 2021;Turner et al., 2022).Adding those effects into Eq.( 19) reflects the effects of formation damage on well injectivity during CO2 storage in aquifers and depleted gas fields.The governing system includes balance equations for energy and mass of salt, methane, and chemical species.In large scale approximation, the hyperbolic system of the quasi-linear equations contains shocks of temperature, saturation, and component concentrations; the corresponding Riemann problems are self-similar (Lake, 1989;Polyanin et al., 2002;Polyanin & Zaitsev, 2003).Accounting for dissipative effects of diffusion and reaction kinetics induces the smoothening of the shocks by travelling wave solutions (Bedrikovetsky, 1993;Poli ︠ a ︡ nin & Dilʹman, 1994).In several cases, the obtained analytical solutions for water displacement allow for exact and asymptotic upscaling (Cheng & Rabinovich, 2021;Rabinovich et al., 2015); numerous references for upscaling relevant to CO2 injection are available from review (Bedrikovetsky & Borazjani, 2022).In addition, the analytical solution can be obtained for a temperature dependency of the J-T coefficient.This will introduce non-linear term   () in Eq. ( 12).For this equation, the characteristics are straight lines.The solution along characteristics is obtained by separation of variables and is given by implicit integral in temperature.
Fig. 9 illustrates how the system parameter ω impacts the validity of Newton's law within the system.As ω decreases, the period during which Newton's law remains valid (tD max) increases, and for the lowest value, there is no such time limit.In general, for every  ≪ 1, there does exist such  and   that for  <   ,   − ( = 0, ) < .Furthermore, Fig. 10 indicates that for injection processes with specific values for A*, V*, and ω, Newton's law holds over a longer observation period when the radius is larger.

CONCLUSIONS
This paper presents an analytical model for a 1D axi-symmetric problem, in which CO2 is injected with constant rate and temperature into a semi-infinite porous medium saturated by gas and water.The model accounts for the Joule-Thomson effect due to "isenthalpic" expansion of the gas and quasi steady-state heat exchange -Newton's law -between the reservoir and the adjacent (over-and under-burden) layers.The exact solution leads to explicit expressions for the temperature profile.The solution shows that the temperature front lags significantly behind the CO2 front due to instant heat exchange between the rock, connate water, and CO2.The exact solution ahead of the front depends on initial temperature and is independent of the injected temperature.In adiabatic reservoir without heat supply from the adjacent layers, the temperature decreases from the well to a minimum value on the heat front, then jumps up on front followed by slow increase up to the initial reservoir temperature.The heat supplied by the adjacent layers has two major impacts on the temperature profile: (1) it decelerates the penetration depth propagation, and (2) it decreases the temperature decline from the well to heat front.Also, the minimum temperature for this case is lower than at the case with no heat gain from the surrounding layers, and its position and value with time remain unchanged.We show that the Newton's heat-exchange law yields tending of the temperature profile to steady-state profile as time tends to infinity; the penetration depth position stabilises.The position of the steadystate temperature strongly depends on dimensionless Joule-Thomson number A*; the dimensionless heat front velocity V* affects less.The results of this study have major implications for CO2 storage in depleted gas fields and can be used as guide to quantify whether the temperature profile in the reservoir falls within the hydrate formation zone, or whether the induced temperature gradient can jeopardize the mechanical integrity of the rock.
Considering the piston-like displacement of water by CO2, two temperature domains are delineated by a shock front.
First, we consider the domain ahead of the heat front, where IC propagates along the characteristic lines.Here the free variable is the dimensionless time,   .
Assuming that,   =  * (  (  ),   ), the total derivative of  * with respect to   is obtained as: Therefore, the domain ahead of heat front is defined by a system of first order ODE as: The system of equations Eq. ( A2) is solved with the following IC: By integrating the first part of the system of equations Eq. (A2), we obtain the trajectory of characteristics as: Further, we solve the second part of the system of equations Eq (A2) to obtain a general solution for temperature in the domain   <   2 −1 2 * as: Here, the () function is defined as: By using the IC Eq. (A3), the constant of integration c1 in Eq. ( A5) is obtained as: By substituting Eq. (A7) into the general solution Eq. (A5), the temperature profile in the domain ahead of the heat front is obtained as: Next, we consider the domain behind of the heat front, where BC propagates along the characteristic lines.Here the free variable is the dimensionless radius,   .Assuming that,   =  ̅ (  ,   (  )), the total derivative of  ̅ with respect to   is obtained as: Therefore, the following system of first order ODE is obtained for the domain behind of the heat front.Finally, the shock emanates from the point (1,0) in the rD -tD plane; which allows for the overall dimensionless solution for the problem defined by Eq. ( 13) with IC Eq. ( 14) and BC Eq. ( 15) to be obtained as:
(i) 1D radial unsteady-state single-phase non-isothermal flow in an infinite homogeneous reservoir; (ii) water, CO2, and rock are incompressible; (iii) the reservoir (with permeability k [L 2 ] and porosity ϕ) contains water with an initial water saturation S=1; piston-like displacement of water by CO2 resides constant saturation Swi of immobile water behind the displacement front; (iv) temperature at outer boundaries of under and overburden shales is equal to initial temperature TI [K]; (v) the temperature of the injected CO2 in the wellbore is constant TJ [K]; (vi) Joule-Thomson coefficient, permeability, porosity, heat conductivity and heat capacity of rock and fluid are constant; (vii) CO2 is injected at a constant rate of qJ [M 3 T -1 ], and pressure at drainage ; (viii)  water evaporation into gaseous phase is neglected; (ix) density and viscosity of fluids are constant; (x) heat exchange between surrounding shales and the reservoir occurs under steady-state conditions.

Fig. 3 .
Fig. 3. Schematic of CO2 injection into a depleted reservoir.The cold part (due to Joule-Thomson cooling effect) is shown in blue.

Fig. 4 .
Fig. 4. Joule-Thomson coefficient of CO2 as a function of pressure for different temperatures (modified from Gao et al. (2021) )

4. 2
Formulating the Criterion for Validity of the Heat Exchange Term.

Fig. 7 .
Fig. 7.A tetrahedron defining the domain of the three dimensionless constants in the system.
1 − (  = 0,   ) against dimensionless time for various radial distances in the reservoir in Fig 8a.We also obtain a generalized area of validity in the (rD, tD) plane in Fig 8b, where rc is the critical radius.

Fig 8 .
Fig 8. Validity of steady state heat exchange (a) validity vs. time for different radii.(b) domain of validity

Fig. 10 .
Fig. 10.Difference between temperature on the outer bounds of shale and initial reservoir temperature, to validate Newton's law, temperature history at the distance from injectors at (a) r = 0.002 rinj, (b) r = 0.01 rinj, (c) r = 0.015 rinj, (d) r = 0.03 rinj Fig. 12.(a) The impact of the heat transfer coefficient, κ, on the temperature profile.(b).The impact of the heat transfer coefficient, κ, on the temperature penetration depth.

Fig. 14 .
Fig. 14.The impacts of the mass injection, permeability and reservoir thickness on the minimum temperature and the steadystate position of the temperature front in the reservoir.The axes are normalized to the parameters of the base case.

Table 1 .
Model parameters used in this study.
Fig.6.Formation of a region with a non-zero water saturation between dry-out and cold fronts.