Numerical Simulation of Hydrogen Diffusion in Cement Sheath of Wells Used for Underground Hydrogen Storage

: The negative environmental impact of carbon emissions from fossil fuels has promoted hydrogen utilization and storage in underground structures. Hydrogen leakage from storage structures through wells is a major concern due to the small hydrogen molecules that diffuse fast in the porous well cement sheath. The second-order parabolic partial differential equation describing the hydrogen diffusion in well cement was solved numerically using the ﬁnite difference method (FDM). The numerical model was veriﬁed with an analytical solution for an ideal case where the matrix and ﬂuid have invariant properties. Sensitivity analyses with the model revealed several possibilities. Based on simulation studies and underlying assumptions such as non-dissolvable hydrogen gas in water present in the cement pore spaces, constant hydrogen diffusion coefﬁcient, cement properties such as porosity and saturation, etc., hydrogen should take about 7.5 days to fully penetrate a 35 cm cement sheath under expected well conditions. The relatively short duration for hydrogen breakthrough in the cement sheath is mainly due to the small molecule size and high hydrogen diffusivity. If the hydrogen reaches a vertical channel behind the casing, a hydrogen leak from the well is soon expected. Also, the simulation result reveals that hydrogen migration along the axial direction of the cement column from a storage reservoir to the top of a 50 m caprock is likely to occur in 500 years. Hydrogen diffusion into cement sheaths increases with increased cement porosity and diffusion coefﬁcient and decreases with water saturation (and increases with hydrogen saturation). Hence, cement with a low water-to-cement ratio to reduce water content and low cement porosity is desirable for completing hydrogen storage wells.


Introduction
Carbon emissions from fossil fuels causing global warming have grown exponentially since the 19th century. As a result, energy companies strive to reduce their carbon footprint by transitioning to large-scale deployment of renewable energy. Although renewables such as solar and wind are intermittent in supply, hydrogen is the most abundant element and a primary sustainable renewable energy source.
The low energy density of hydrogen implies that it requires more volumetric storage capacity than natural gas to produce the same energy [1,2]. Consequently, large-scale subsurface porous media, such as salt caverns, aquifers, and depleted oil and gas reservoirs, are more practical options for underground hydrogen storage (UHS) [3].
Depleted oil and gas reservoirs are the most sorted candidates for UHS; they have rich databases to support Underground Storage (UGS) site activities and have been studied extensively by geoscientists and reservoir engineers. In addition, studies show that the geological trap mechanisms of most depleted oil and gas reservoirs demonstrate structural integrity and safety, as they have kept their hydrocarbons for millions of years without leaks. Finally, the economics of injecting and producing hydrogen from depleted oil and gas diffusion decreases with increasing relative humidity. Teodoriu et al. [31] studied the effects of salt concentration on class G cement for wells used in underground gas storage sites. The study analyzed the thickening time, compressive strength, elastic modulus, and set cement permeability of class G cement at different salt concentrations. The researchers were curious about how salts affected the performance and integrity of cement used in underground gas storage wells by varying the salt concentration.
The numerical simulation of hydrogen migration into cement sheaths is an essential tool for researchers and engineers to understand the complex processes involved in the degradation of oil and gas well cementing in UHS. It provides a powerful approach to modeling the physical and chemical phenomena involved in hydrogen migration and predicting the behavior of cement sheaths under different scenarios. Ghaedi et al. [32] used a series of one-and three-dimensional numerical models to investigate hydrogen diffusion through the caprock during geological storage. In addition, a semi-analytical solution was proposed to predict the cumulative hydrogen loss over time.
Unfortunately, there are no or limited numerical simulations of hydrogen diffusion through the cement sheathes of oil and gas wells used in UHS. Hence, the primary goal of this study was to provide numerical simulations of hydrogen diffusion through the cement sheath of a well used in hydrogen storage by simulating Fick's second law of diffusion for fixed-concentration boundary conditions. In addition, a modified effective diffusion coefficient of hydrogen was introduced into Fick's second law of diffusion to account for the presence of other substances in the cement pores. Simulation models were then verified using an analytical model. Sensitivity analysis was performed to understand better the complex processes involved in the degradation of the cement sheath with invariant matrix and fluid properties by hydrogen migration. Finally, this work simulated the hydrogen migration time along the axial direction of a cement column with varying matrix and fluid properties above the hydrogen storage zone.

Mathematical Model
Fick's second law of diffusion is a generalization of Fick's first law of diffusion, which states that the flux of a substance across a plane is proportional to the concentration gradient in that plane [33][34][35][36][37].
Fick's second law of diffusion is a second-order parabolic partial differential equation describing a substance's diffusion in a medium. It considers the substance's time-dependent (non-steady state) diffusion in the medium. Based on this study, Fick's second law of diffusion will describe the diffusion of hydrogen (substance) through a cement sheath (medium).
The solution to Fick's second law of diffusion depends on the boundary and initial conditions. For fixed-concentration boundary conditions, the concentration of hydrogen is fixed at one end of a one-dimensional medium of length L, and there is initially no hydrogen in the medium.
Fick's second law of diffusion is given by: Initial Condition: Boundary Conditions: where C(x, t) is the concentration of the substance (hydrogen), t is time, x or L is the distance (cement thickness), S w is the saturation of water in cement pores, (1 − S w ) is the saturation of gas (hydrogen) in the cement pore, Ø is the porosity of the cement sheath, and D e f f is the effective diffusion coefficient of hydrogen in the presence of other substances in the cement pores, while D o is the absolute diffusion coefficient of hydrogen in the absence of other substances in the cement pores. The diffusion coefficient will differ depending on subsurface variables such as concentration, temperature, and pressure [35]. Equation (1), subject to Neumann boundary conditions, can be solved analytically using Laplace transformations [38][39][40]. The analytical solution is presented below: The series converges fast, even at n ≤ 10.
Application of this solution is limited to homogeneous systems where both the porous matrix and the fluid are invariant in space. But in reality, matrix and fluid properties can vary spatially. For example, the porosity of well cement varies with depth, and the water saturation in the pore space also changes; in such conditions, a numerical model is highly desirable. The Finite Difference Method (FDM) was used to solve the equations numerically, and solutions were implemented with Python. The numerical formulation with second-order accurate approximation is shown in Appendix A. Three assumptions made in the numerical model are: (1) hydrogen is non-dissolvable in the water present in cement pores; (2) the diffusion coefficient of hydrogen remains constant in the UHS; and (3) the cement properties, such as porosity and saturation, could be constant or variable depending on the problem setup.

Model Validation
This section aims to investigate the correctness of the numerical model. We verify the numerical model presented in Appendix A with the analytical solution given by Equation (6). Table 1 presents input data to evaluate the hydrogen diffusion into the cement sheath of a well in a UHS site. The data were obtained from literature on well completion and SHASTA/DOE 2022 report. Initial normalized concentration of hydrogen, C 0 1 Figure 1 compares the numerical results and the analytical solution after 2.89 days of hydrogen diffusion, demonstrating a good match between the two solutions. The solid blue curve represents the numerical solutions, while the red dotted curve represents the analytical solution. However, one challenge with the study is the absence of experimental studies, such as Scanning Electron Microscope/Energy Dispersive Spectroscopy (SEM-EDS) mapping and SEM-EDS point analysis for investigating the depth of hydrogen penetration in class G or H well cement under UHS conditions for further validation of the numerical model.

Numerical Experiments and Analysis of Results
While the absence of experimental data may limit the direct validation of our numerical studies, we employed alternative approaches, such as sensitivity analysis. Sensitivity analyses were conducted to account for uncertainties in some parameter values presented in Table 1, such as coefficient of diffusion in the cement, cement porosity, gas saturation, and cement thickness.
In Figures 2-4, simulation time ranging from 0.12 days to 69.44 days was run to expose the cement sheath to hydrogen diffusion until a breakthrough occurred. Hydrogen took about 7.5 days to penetrate the cement sheath fully. The minimum time for full penetration is when the hydrogen concentration at the end of the cement thickness reaches 10% of the initial normalized hydrogen concentration. Figures 2-4 also show that it will take about 69 days for 90% of the initial hydrogen concentration to penetrate the cement sheath fully.

Numerical Experiments and Analysis of Results
While the absence of experimental data may limit the direct validation of our numerical studies, we employed alternative approaches, such as sensitivity analysis. Sensitivity analyses were conducted to account for uncertainties in some parameter values presented in Table 1, such as coefficient of diffusion in the cement, cement porosity, gas saturation, and cement thickness.
In Figures 2-4, simulation time ranging from 0.12 days to 69.44 days was run to expose the cement sheath to hydrogen diffusion until a breakthrough occurred. Hydrogen took about 7.5 days to penetrate the cement sheath fully. The minimum time for full penetration is when the hydrogen concentration at the end of the cement thickness reaches 10% of the initial normalized hydrogen concentration. Figures 2-4 also show that it will take about 69 days for 90% of the initial hydrogen concentration to penetrate the cement sheath fully.

Numerical Experiments and Analysis of Results
While the absence of experimental data may limit the direct validation of our numerical studies, we employed alternative approaches, such as sensitivity analysis. Sensitivity analyses were conducted to account for uncertainties in some parameter values presented in Table 1, such as coefficient of diffusion in the cement, cement porosity, gas saturation, and cement thickness.
In Figures 2-4, simulation time ranging from 0.12 days to 69.44 days was run to expose the cement sheath to hydrogen diffusion until a breakthrough occurred. Hydrogen took about 7.5 days to penetrate the cement sheath fully. The minimum time for full penetration is when the hydrogen concentration at the end of the cement thickness reaches 10% of the initial normalized hydrogen concentration. Figures 2-4 also show that it will take about 69 days for 90% of the initial hydrogen concentration to penetrate the cement sheath fully.  The relatively short duration for hydrogen breakthrough in the cement sheath is mainly due to the small molecule size and high hydrogen diffusivity. Therefore, the well cement integrity is paramount to ensure any successful UHS project.    The relatively short duration for hydrogen breakthrough in the cement sheath is mainly due to the small molecule size and high hydrogen diffusivity. Therefore, the well cement integrity is paramount to ensure any successful UHS project.    The relatively short duration for hydrogen breakthrough in the cement sheath is mainly due to the small molecule size and high hydrogen diffusivity. Therefore, the well cement integrity is paramount to ensure any successful UHS project.  Figure 5 shows the sensitivity of the coefficient of diffusion of hydrogen. A hydrogen diffusion coefficient ranging from 0.019 cm 2 /s to 0.001945 cm 2 /s was selected to account for uncertainties in the actual hydrogen diffusion coefficient in cement used in UHS. Based on the SHASTA/DOE 2022 report [6], the diffusion coefficient decreased with curing time and increased with water-cement ratio. For example, for a water-cement ratio of 0.5, the diffusion coefficient was 0.019 cm 2 /s after 28 days of curing. We observed that an increase in the hydrogen diffusion coefficient would cause a corresponding increase in the hydrogen diffusion rate in the cement sheath. Conversely, decreasing the diffusion coefficient will result in a slower diffusion rate. In Figure 5, when the diffusion coefficient is in the magnitude of 10 −2 cm 2 /s, hydrogen will fully penetrate the cement sheath within one to two days, but when the diffusion coefficient is in the magnitude of 10 −3 cm 2 /s it takes a few months for hydrogen will fully penetrate the cement sheath. It is important to note that the diffusion coefficient could differ depending on subsurface variables such as concentration, temperature, and pressure.  Figure 5 shows the sensitivity of the coefficient of diffusion of hydrogen. A hydrogen diffusion coefficient ranging from 0.019 cm 2 /s to 0.001945 cm 2 /s was selected to account for uncertainties in the actual hydrogen diffusion coefficient in cement used in UHS. Based on the SHASTA/DOE 2022 report [6], the diffusion coefficient decreased with curing time and increased with water-cement ratio. For example, for a water-cement ratio of 0.5, the diffusion coefficient was 0.019 cm 2 /s after 28 days of curing. We observed that an increase in the hydrogen diffusion coefficient would cause a corresponding increase in the hydrogen diffusion rate in the cement sheath. Conversely, decreasing the diffusion coefficient will result in a slower diffusion rate. In Figure 5, when the diffusion coefficient is in the magnitude of 10 −2 cm 2 /s, hydrogen will fully penetrate the cement sheath within one to two days, but when the diffusion coefficient is in the magnitude of 10 −3 cm 2 /s it takes a few months for hydrogen will fully penetrate the cement sheath. It is important to note that the diffusion coefficient could differ depending on subsurface variables such as concentration, temperature, and pressure.  Figure 6 shows the sensitivity of hydrogen saturation on the cement sheath. The result shows that when hydrogen is exposed to the cement sheath for a given period, the hydrogen diffusion rate in the cement sheath is higher when the water saturation in the cement pores is lower (higher hydrogen saturation) and vice versa. Figure 6 also shows that when the water saturation in the cement pore is 40% (i.e., 60% hydrogen saturation), the normalized hydrogen concentration in the fully penetrated cement sheath is 0.6. Conversely, when the water saturation in the cement pore is 90% (i.e., 10% hydrogen saturation), the normalized hydrogen concentration in the penetrated cement sheath is about 0.02. Figure 7 shows the sensitivity of the cement sheath porosity. When hydrogen is exposed to the cement sheath for a given period, the hydrogen concentration in the penetrated cement sheath increases with increased cement sheath porosities. Hence, cement with a low water-to-cement ratio, high compressive strength, and low cement porosity is desirable for wells used in UHS.  Figure 6 shows the sensitivity of hydrogen saturation on the cement sheath. The result shows that when hydrogen is exposed to the cement sheath for a given period, the hydrogen diffusion rate in the cement sheath is higher when the water saturation in the cement pores is lower (higher hydrogen saturation) and vice versa. Figure 6 also shows that when the water saturation in the cement pore is 40% (i.e., 60% hydrogen saturation), the normalized hydrogen concentration in the fully penetrated cement sheath is 0.6. Conversely, when the water saturation in the cement pore is 90% (i.e., 10% hydrogen saturation), the normalized hydrogen concentration in the penetrated cement sheath is about 0.02. Figure 7 shows the sensitivity of the cement sheath porosity. When hydrogen is exposed to the cement sheath for a given period, the hydrogen concentration in the penetrated cement sheath increases with increased cement sheath porosities. Hence, cement with a low water-to-cement ratio, high compressive strength, and low cement porosity is desirable for wells used in UHS.  One advantage of the numerical simulation approach to hydrogen diffusion through the cement sheath over the analytical solution is the ability to simulate varying properties of cement, such as distributed porosity and water saturation in the pore spaces. For example, when investigating the time taken for hydrogen migration in the axial direction of the cement column above the storage zone (50 m of cement along the caprock thickness and 100 m above the caprock), if the cement column's pore spaces have a varied distribution of cement porosities and water saturations and no crack on the cement sheath, as shown in Figure 8 below, then the time taken to penetrate the 50 m cement column across the caprock is estimated to be about 500 years.  One advantage of the numerical simulation approach to hydrogen diffusion through the cement sheath over the analytical solution is the ability to simulate varying properties of cement, such as distributed porosity and water saturation in the pore spaces. For example, when investigating the time taken for hydrogen migration in the axial direction of the cement column above the storage zone (50 m of cement along the caprock thickness and 100 m above the caprock), if the cement column's pore spaces have a varied distribution of cement porosities and water saturations and no crack on the cement sheath, as shown in Figure 8 below, then the time taken to penetrate the 50 m cement column across the caprock is estimated to be about 500 years. One advantage of the numerical simulation approach to hydrogen diffusion through the cement sheath over the analytical solution is the ability to simulate varying properties of cement, such as distributed porosity and water saturation in the pore spaces. For example, when investigating the time taken for hydrogen migration in the axial direction of the cement column above the storage zone (50 m of cement along the caprock thickness and 100 m above the caprock), if the cement column's pore spaces have a varied distribution of cement porosities and water saturations and no crack on the cement sheath, as shown in Figure 8 below, then the time taken to penetrate the 50 m cement column across the caprock is estimated to be about 500 years.

Conclusions
In this study, the second-order parabolic partial differential equation describing hydrogen diffusion in the cement sheathes of wells used in UHS was solved numerically using the finite difference method (FDM) and implemented with Python. The numerical solution was verified with an analytical solution for an ideal case where matrix and fluid have invariant properties.
The significance of the numerical simulation approach over the analytical method is the ability to simulate hydrogen diffusion in cement with varying properties, such as geometry, distributed porosity, and water saturation in the pore space. This numerical approach will prove useful in expanding our understanding of how wells used for UHS can best be selected, monitored, and managed. This study is limited by the absence of experimental studies, such as Scanning Electron Microscope/Energy Dispersive Spectroscopy (SEM-EDS ) mapping and SEM-EDS point analysis, for investigating the depth of hydrogen penetration in class G or H well cement under UHS conditions and as a result the numerical solution was not validated experimentally.
The following conclusions are drawn based on the sensitivity analysis conducted in this study: 1. Hydrogen should take about 7.5 days to fully penetrate a 35 cm cement sheath under the condition in Table 1. The relatively short duration for hydrogen breakthrough in the cement sheath is mainly due to the small molecule size and high hydrogen diffusivity. If the hydrogen reaches a vertical channel behind the casing, a hydrogen leak from the well is soon expected. Therefore, the well cement integrity is paramount to ensuring any successful UHS project. 2. Hydrogen diffusion into cement sheaths increases with increased cement porosity and diffusion coefficient and decreases with water saturation (increases with hydrogen saturation). Hence, cement with a low water-to-cement ratio to reduce water content and low cement porosity is desirable for wells used in UHS. 3. The simulation result implies that the time taken for hydrogen migration along the axial direction of a 50m cement column across the caprock and above the storage zone is estimated to be about 500 years. Hence, provided the cement sheath is crack-

Conclusions
In this study, the second-order parabolic partial differential equation describing hydrogen diffusion in the cement sheathes of wells used in UHS was solved numerically using the finite difference method (FDM) and implemented with Python. The numerical solution was verified with an analytical solution for an ideal case where matrix and fluid have invariant properties.
The significance of the numerical simulation approach over the analytical method is the ability to simulate hydrogen diffusion in cement with varying properties, such as geometry, distributed porosity, and water saturation in the pore space. This numerical approach will prove useful in expanding our understanding of how wells used for UHS can best be selected, monitored, and managed. This study is limited by the absence of experimental studies, such as Scanning Electron Microscope/Energy Dispersive Spectroscopy (SEM-EDS ) mapping and SEM-EDS point analysis, for investigating the depth of hydrogen penetration in class G or H well cement under UHS conditions and as a result the numerical solution was not validated experimentally.
The following conclusions are drawn based on the sensitivity analysis conducted in this study:

1.
Hydrogen should take about 7.5 days to fully penetrate a 35 cm cement sheath under the condition in Table 1. The relatively short duration for hydrogen breakthrough in the cement sheath is mainly due to the small molecule size and high hydrogen diffusivity. If the hydrogen reaches a vertical channel behind the casing, a hydrogen leak from the well is soon expected. Therefore, the well cement integrity is paramount to ensuring any successful UHS project.

2.
Hydrogen diffusion into cement sheaths increases with increased cement porosity and diffusion coefficient and decreases with water saturation (increases with hydrogen saturation). Hence, cement with a low water-to-cement ratio to reduce water content and low cement porosity is desirable for wells used in UHS.

3.
The simulation result implies that the time taken for hydrogen migration along the axial direction of a 50m cement column across the caprock and above the storage zone is estimated to be about 500 years. Hence, provided the cement sheath is crack-free, hydrogen leaks to the surrounding formations may likely happen along the traverse direction of the cement sheath during the operating life of the UHS sites. Acknowledgments: I thank Boyun Guo and Yin Feng for their invaluable guidance and support throughout this publication.

Conflicts of Interest:
The authors declare no conflict of interest.
The implicit FDM are unconditionally stable, converging for all ∆x and ∆t. Also, the above algorithm for Fick's second law of diffusion is implemented in Python.