Analysis of CO 2 Migration in Horizontal Saline Aquifers during Carbon Capture and Storage Process

: The storage of CO 2 has become an important worldwide problem, considering that an excess of CO 2 in the Earth’s atmosphere causes dramatic changes in its climate. One possible solution is to remove the excess of CO 2 from the atmosphere, capture it in the process of creation, and store it safely, negating the possibility of its return into the atmosphere. This is the process of Carbon Capture and Storage (CCS). In the following paper, the authors investigate horizontal saline aquifers and their ability to store CO 2 . The authors’ application of sensitivity analysis on horizontal migrations uncovered that CO 2 permeability and aquifer porosity have a considerable impact on horizontal migrations. During the migration process, CO 2 can reach tens of kilometers from its injection point. By introducing effective CO 2 density to the conduction velocity term, the authors showcase that the convection-diffusion equation for compressible ﬂuids can be replaced with the equation for incompressible ﬂuids. The buoyancy factor in convective velocity is as density dependent as in conduction velocity. By means of introducing an effective density to the aforementioned term, the process of transport via variable convective velocity can be substituted for a process which is effective, constant, and not density dependent.


Introduction
Despite being a trace element, carbon dioxide is an important component of the Earth's atmosphere [1]. It is a key factor of the planet's greenhouse effect [1,2]. The issue with atmospheric CO 2 is that today's concentrations are 150% of its value, which is more than half of that a century ago, and they are still rising. By 2005, the CO 2 concentration increased from around 280 ppm of the 19th century to 380 ppm [2], and its concentration in atmospheric gases is currently around 420 ppm [3]. Today's CO 2 concentrations are the highest they have been in 14 million years [4] and these recorded values are directly attributed to the acts of human beings. Two main factors arise as probable causes, namely the reduction in forest areas and the combustion of fossil fuels. In 2015, the emission of CO 2 into the atmosphere caused by human action was estimated at around 37 Gt/year [5]. The atmospheric concentration of CO 2 is in dynamic balance with the CO 2 dissolved in seas and oceans, and increased levels of CO 2 shift the balance at the junction of hydrosphere, thereby dissolving CO 2 into oceans [6] and reducing its pH value by creating carbonic acid.
The created CO 2 can be removed from the atmosphere naturally via carbon cycles, mostly during the process of plant photosynthesis. The Industrial Age created greater CO 2 production, establishing new and increased values. Natural carbon cycles are not able to reduce CO 2 , and if action is not taken to reduce CO 2 emissions by human and remove the excess from the Earth's atmosphere, future projections do not look promising. Such action of reducing CO 2 is present in the Carbon Capture and Storage (CSS) process [2].
In the pre-combustion route, fossil fuels react with oxygen (O 2 , usually from the air) to produce syngas [13], which is a mixture of CO and hydrogen (H 2 ). Afterwards, syngas is usually purified to remove impurities [14], and then is subjected to a reaction with water vapor, giving CO 2 and H 2. As a result of the precombustion route, CO 2 and hydrogen fuel are the products, where CO 2 can be separated through various processes for storage or utilization [15,16]. The pre-combustion route can have lower costs compared to other methods [17], but its implementation in current facilities have limits. Oxyfuel combustion implies the combustion of fossil fuels with pure oxygen, O 2 . This method is the most promising one but requires high-cost pure oxygen. Cryogenic air separation is one of the most utilized methods [18] but requires a large amount of energy for cryogenic air separation of O 2 . The post-combustion route involves the capture and separation of CO 2 from combusted gas [19,20]. Gas created in combustion processes needs to be cooled, removed from dust particles, and purified [21].
Post-combusted technologies are mostly used, since these processes appeared in the beginning of CCS processes and are more developed [22,23]. Capture technologies incorporate various techniques for CO 2 separation, such as absorption, separation with membranes, and adsorption [24]. These techniques are based on physical, chemical, and even biological systems.

Carbon Storage
Carbon storage, or carbon sequestration, is the process of storing CO 2 , long-term or permanently, using various methods, such as dissolving it in the ocean or using geological storage [25]. Geological storage has the highest number of long-term benefits [26], and it comes in the form of two options: depleted oil-gas reservoirs and saline aquifers. Furthermore, CO 2 can be injected into oil reservoirs that are still in use as a form of enhanced oil recovery [27][28][29].
The advantage of exhausted oil and gas reservoirs lies, first of all, in their ability to unequivocally store CO 2 ; since reservoirs are capable of holding oil and gas for millions of years, CO 2 can remain trapped for long periods of time. The only way for CO 2 to escape is through the drilling point, which should be sealed properly. The properties of such reservoirs are known due to hydrocarbon extraction, when the amount of equipment for measuring and tracking, installed on site, is the highest [30]. Regardless, the disadvantages of such reservoirs lie in their storage capacities, especially those of gas reservoirs. If the densities of natural gas, oil, and underground CO 2 in average reservoir conditions are 150, 800, and 700 kg/m 3 , respectively, then 1 m 3 of the gas can be replaced with 0.59 m 3 of CO 2 and 1m 3 of the oil can be replaced with 2.75 m 3 of CO 2 [31]. However, 1 kg of natural gas produces 2.75 kg of CO 2 , while 1 kg of oil produces 3.14 kg of CO 2 [30]. Burning oil, therefore, produces more CO 2 than can be stored at the site where oil is drilled.
Storing CO 2 in deep saline aquifers is arguably the best option now [31][32][33][34][35][36][37]. In contrast to the reservoirs depleted from gas and oil, aquifers have larger storage capacities and exist in great abundance worldwide. These formations can be in the form of layers of permeable rocks, such as sediments, surrounded at the top with less permeable caprocks. The permeable layer is filled with saline water, and when CO 2 is pumped into saline aquifers, it replaces the water and tends to diffuse into porous rocks, where it is ultimately captured [38]. Considering that upward diffusion needs to be prevented, the top layer of a saline aquifer caprock [39] must be from a material which has a low permeability, such as clay [34]. The diffusion of CO 2 in such a condition is, due to porous material, enabled horizontally but suppressed upwardly. The process is called structural trapping and it increases the storage capacities of saline aquifer reservoirs.

Carbon Escape from Structural Storage
The buoyancy force tends to steer CO 2 into vertical directions, thereby leaking CO 2 through the cracks in caprock formations that lie above the aquifers [38,39]. Makhnenko and coauthors [35] concluded that caprock microstructures strongly influence the pressure and permeability of CO 2 and its ability to escape through caprock. To plan geological storage projects, the International Energy Agency [40] suggests that for each formation, the evolution of its properties in interaction with CO 2 should be studied [41]. The effects of interaction can be quite opposite. Some cases show the sealing of pores and decreasing permeability [42,43], while in other cases studies showed an increase in permeability because of the interaction with CO 2 [44,45]. Pumped CO 2 is often not in a supercritical state, but in a liquid state. This implies that the temperature of pumped CO 2 is far less than the temperature in the injection point, leading to non-equilibrium state in which CO 2 can flows downwards along the injection well [46]. This can expand the spreading of CO 2 in lateral directions and increase the chance of its escape if it reaches the caprock endings. For long-scale injections, the constant pumping of cooler CO 2 can cause cooling of the caprock, leading to changes in the physical parameters of the caprock [47][48][49][50][51][52][53].
Considering that CO 2 injection is usually performed in sedimentary basins, oil and gas extractions may lead to drilling sites and exploration wells in abandon basins, causing CO 2 to escape [54]. The oil exploration process has generated a vast number of explorations wells, which exist for the purpose of finding oil in sedimentary basins. Hence, dry wells without seals and precautions are large in number [55]. Depending on the location, well densities can be at the order of several wells per square kilometer [49]. If CO 2 can migrate laterally in a radius of, for example 5 km, there could be hundreds of existing wells. Moreover, wells with concrete sealing that are centuries old have been degraded, presenting openings through which CO 2 can easily escape the aquifer.
It must be ensured that residual trapping occurs during the storage process, i.e., the entrapment of CO 2 between an impermeable layer and saline water or structural trapping [50]. Physically trapped CO 2 leads to the second stage of CO 2 trapping called capillary trapping in which bulk CO 2 plume is divided into CO 2 blobs, disconnected due to the depression in the caprock [51,52]. The dissolution of CO 2 in aquifer brine follows the process called solubility trapping [56][57][58][59], followed by long-term CO 2 mineralization [60,61] or mineral trapping. Mineralization is crucial for trapping CO 2 because it ensures longterm immobilization, thereby reducing the probability of CO 2 escaping and returning into the atmosphere [62]. Still, CO 2 dissolution is a more dominant storage mechanism than mineralization [63]. Solubility trapping can show significant results one hundred years from the point of injection [64], whereas mineralization takes far longer [65,66].
Both oil reservoirs and saline aquifers have depth distributions of 700-800 m below the Earth's ground level [30,67,68]. Saline aquifers with depths below 800 m have temperatures of 40-50 • C and pressures of up to 8 MPa [30,67]. In depths that reach below 800 m, CO 2 can be in a supercritical state, with a temperature above 31 • C and the critical pressure of 7.3773 MPa. In such a supercritical state, CO 2 is not liquid nor gas, but possesses the properties of both. It's density is equivalent to that of a liquid, but it exhibits diffusion and viscosity as if it were a gas, implicating that it dissolves in water like a liquid, but diffuses through solids like gas. Moreover, the aquifers' storage capacity depends on their parameters, and the density of CO 2 in a supercritical state can range from 150 to 800 kg/m 3 and higher [68].
In the literature, papers devoted to the problem of CO 2 storage and its migration in saline aquifers deal with plume formation in the injection period [40, [68][69][70]. Lateral brine migrations need to be concerned in case of the potential pressure driven by brine leakage through localized pathways, such as unsealed drilling holes, faults, and less component caprocks [71]. The horizontal/lateral diffusion is enabled in such conditions and the trap is referred to as a structural trap. CO 2 can leak through structural traps by escaping through the cracks or pores of caprocks with low permeabilities or through artificial pathways such as surrounding wells [2]. When it comes to saline aquifers, CO 2 is injected into formations that are at the top of saline water reservoirs. When equilibrium is established at the interface between water and CO 2 's supercritical phase of equilibrium, CO 2 dissolves into the aquifer. The process is influenced by the aquifer's parameters, the main ones being salt concentration, temperature, and pressure, and will differ for every location. Much effort has been dedicated to modeling real-time conditions in saline aquifers and determining the values of parameters which influence the dissolubility of CO 2 and its migration pathways [72][73][74][75]. Recent approaches have entailed injecting CO 2 foam so as to reduce the mobility of CO 2 and trap it at a more successful rate [76].
The diffusion rate of CO 2 is an important parameter for estimating its dissolubility and reports from experiments are mainly devoted to its diffusivity in pure water, which has been determined to range from 1.15 to 8.2 × 10 9 m 2 /s for temperature ranges from 10 to 90 • C [56]. Azin [77] recreated such conditions in a laboratory in order to determine CO 2 's diffusion coefficient in saline aquifers with pressures at 6-7 MPa, finding that the diffusion coefficient ranges from 3.5 × 10 −9 to 6 × 10 −9 m 2 /s. Fick's second law has also been utilized to model CO 2 inside geological storages, along with additional terms which evaluate the mass balance of the gas phase [77].

Research Goal
We assume that some amount of CO 2 is injected into a structural trap. Supercritical CO 2 , surrounded by denser water, forms a bulb or structural brine in the caprock celling during the first phase and starts to diffuse laterally in the second phase. Relying on the injection point and the amount of CO 2 inserted into the aquifer during the first phase, the buoyancy force and difference in pressure guides the flow of the fluid which has gone through two phases. CO 2 with a lesser density tends to replace water and form an equilibrium structure at the top of the aquifer, beneath the caprock. This formation then starts the process of lateral migration. The aim of this article is to examine the ways in which CO 2 can escape due to lateral migration pathways inside aquifers, where residual trapping is not fully possible. Worst-case scenario estimation is performed in the case of lateral migrations.
In the migration process, CO 2 can reach abundant wells which are not properly sealed or not sealed at al. Furthermore, CO 2 can reach some natural gaps in the caprock structure or some porosity layers which enable upwards migrations. A further step is the exploration of the possibility of upward migration and leakage in the worst-case scenario. Additional research is performed in this manuscript to estimate the worst-case scenario breach of the CO 2 from saline aquifers and leaking into the atmosphere.
For the purposes of this research, mathematical methods of CO 2 transport in lateral and vertical directions were used based on the transport equations-diffusion and convection. In this manuscript, the authors also explore the possibility of finding the simplification of mathematical models without of losing generality and accuracy in the description of transport processes of CO 2 .

Materials and Methods
To model this flow of CO 2 and determine the range of lateral migration in saline aquifers, Nordbotten [78] modeled the sharp interface between CO 2 and brine (saline water). No mixture, as a result of diffusing CO 2 , is of interest since this process of migration can be neglected compared to fluids like the flow of supercritical CO 2 . The diffusion process has its role later.
Vertically averaged pressure which governs the flow inside the aquifer, wherein CO 2 is injected with the rate of Q and at the depth of H below the caprock, can be mathematically modeled as the pressure gradient [78] −kH∇(λ∇p) = Q · δ(r − r'), where k presents permeability, λ is the total mobility, δ denotes the Dirac delta, r is the coordinate of fluid flow, and r' is the injection coordinate. The fluid's mobility is defined as the ratio of the relative permeability and fluid viscosity µ. Total mobility λ is defined as the mobility λ c of CO 2 and λ w of water where h is the thickness of the CO 2 phase. The equation that describes the evolution of the sharp interface [78][79][80] or the plume is given in Appendix A-Equation (A1). Nordbotten [78] derived the thickness of h from the CO 2 plume as a function of the radial distance and time: wherein ϕ is the porosity of the aquifer and V(t) is the total volume of CO 2 inserted up to a time of t. This can be determined from the injection rate V(t) = Q(t)dt. For a constant injection rate, volume is the linear function of time V(t) = Qt. In an arbitrary injection period, the total amount of volume injected is V tot , and the maximum radial distance that CO 2 reaches during the injection period can be determined from Equation (1) by setting the thickness of the plume to zero r max (t) = λ c V tot ϕπλ w H . The parameters needed to determine radial flow of the CO 2 are location dependent. The depth of the aquifer, pressure, temperature, porosity and other parameters influence the radial flow. The pressure gradient can be taken to be around 10.5 MPa/km, while the temperature gradient can vary from 25 • C/km up to [70] and for depths around 1000 m the viscosity of the CO 2 and brine is in range of µ c = (0.0611 − 0.0395)mPa · s and µ w = (0.195 − 0.644)mPa · s, respectively. Depending on the temperature of the basins in which the aquifers are found, the density of the CO 2 phase and water can significantly vary. In deep formations, the CO 2 density can be from 733 kg/m 3 in cold basins down to 479 kg/m 3 in warm basins. The water density can vary from 1202 kg/m 3 down to 94 kg/m 3 [70].
The second stage of CO 2 migration ensues after the forming of equilibrium, wherein the shape and radial distance from the injection site are derived from the properties of the aquifer. Assuming that CO 2 is permanently captured, its dissolving into water becomes predominant due to diffusion and mineralization [81]. While these processes are present from the start of the injection, they are not as present compared to flow processes, which shape the formation of the brine. The process of permanently capturing CO 2 is slowmoving and requires periods of geological time, measured in millions of years [81]. If CO 2 reaches the possibility of upward migration during the storage process or the formation of brine, leakage can occur [82]. Therefore, storing CO 2 is not contingent only on volumetric capacities, but on the temporal presence of CO 2 from the point of injection. Leakage can be caused by various processes, mainly by diffusion and viscous flow, which are led by hydrostatic pressure and bouncy force. The 10 m thickness of CO 2 raises the hydrostatic pressure to around 0.03 MPa [81]. While the process takes time, CO 2 can reach shallow, drinking-water basins via upward motion and, ultimately, surface from the ground, thereby leaking into the atmosphere. Carbon sequestration is a long-term process and CO 2 storage should, hence, be planned long ahead. This article attempts to explore the limitations of carbon sequestration for long time periods by addressing mathematical models governed by physical processes for cases where the breach of the CO 2 occurs.
Aside from migrating horizontally, CO 2 can migrate upward, thereby finding pathways into the atmosphere. The most common pathways are cracks inside caprocks and abundant wells, especially if they are not properly sealed [82]. The permeability of seals influence the pressure buildup and brine mobility inside of the storage formation. If the seals have high permeabilities, they can allow leakage vertically upward [82]. Moreover, even the probability of CO 2 escaping through the cement-rock interface in a sealed well is great [83][84][85]; due to the degradation of cement over time, the permeability of this junction can increase in magnitude. The escape rate of CO 2 depends on the type of ground material; for sandstone, porosity can be around 20% and permeability can be at 100 mD, while the Sustainability 2023, 15, 8912 6 of 17 parameters for limestone are 8.5% and 1.5 mD, respectively [86]. Cao, in a paper from 2016 [85], recreated these conditions in a cement-rock interface and found a permeability increase of 3 times for sandstone and up to 24 times for limestone for experiment periods lasting up to 100 h. Furthermore, disruptions in the junction can be even more exacerbated after longer periods of time. Leakage via abundant and improperly sealed wells does not only entail the return of CO 2 into the atmosphere [87], but also the contamination of clean water aqueducts.
To simulate upward migration, the geometry of the problem needs to be defined. Hence, one must consider the injection point as defined in the coordinate system in Figure 1, wherein the origin is the point of injection. An improperly sealed well is at a distance CO 2 can reach, along the r w -axis, and is also geometrically defined in Figure 1. The origin of another coordinate system is placed at a point where CO 2 leaks through a pathway in caprock sealing or rather where the well is established. This coordinate system is used to describe the upwards motion of CO 2 . The ground surface has a vertical coordinate equal to the depth of the aquifer, i.e., H. and abundant wells, especially if they are not properly sealed [82]. The permeability of seals influence the pressure buildup and brine mobility inside of the storage formation. If the seals have high permeabilities, they can allow leakage vertically upward [82]. Moreover, even the probability of CO2 escaping through the cement-rock interface in a sealed well is great [83][84][85]; due to the degradation of cement over time, the permeability of this junction can increase in magnitude. The escape rate of CO2 depends on the type of ground material; for sandstone, porosity can be around 20% and permeability can be at 100 mD, while the parameters for limestone are 8.5% and 1.5 mD, respectively [86]. Cao, in a paper from 2016 [85], recreated these conditions in a cement-rock interface and found a permeability increase of 3 times for sandstone and up to 24 times for limestone for experiment periods lasting up to 100 h. Furthermore, disruptions in the junction can be even more exacerbated after longer periods of time. Leakage via abundant and improperly sealed wells does not only entail the return of CO2 into the atmosphere [87], but also the contamination of clean water aqueducts.
To simulate upward migration, the geometry of the problem needs to be defined. Hence, one must consider the injection point as defined in the coordinate system in Figure  1, wherein the origin is the point of injection. An improperly sealed well is at a distance CO2 can reach, along the rw-axis, and is also geometrically defined in Figure 1. The origin of another coordinate system is placed at a point where CO2 leaks through a pathway in caprock sealing or rather where the well is established. This coordinate system is used to describe the upwards motion of CO2. The ground surface has a vertical coordinate equal to the depth of the aquifer, i.e., H.  The equation that governs the migration of CO 2 is the convection-diffusion equation for compressible fluids, which depicts migration through a crack in the caprock and its junction with cement, sandstone, or limestone: Herein, ρ is the CO 2 density, D is the diffusion coefficient of CO 2 , ϕ is the porosity, and v z is the superficial velocity equal to [88] . Via inserting this velocity into Equation (2), the following is obtained: Equation (3) is the convection-diffusion equation in which density depends on the depth. To solve Equation (3), the finite difference method is employed in this paper. The diffusion-convection equation can be transformed in the discrete form as follows: The backward scheme was used for the first derivate, and the central scheme was used for the second derivate , considering ∆z and ∆t as small quantities. Index i denotes the z coordinate, while index j denotes the time (t). For this equation in the discrete form, the parameters of discretizing space and time must be set, and their values depend on the particular case and calculation. To ensure convergence to the correct values, the condition for numerical solution stability regarding diffusion is [89]: Because of this condition, time domain discretization and space domain discretization are not independent. By setting any of the two parameters ∆z or ∆t, the other is not arbitrary but should follow Equation (5). For example, if the diffusion coefficient is D ≈ 10 −10 m 2 /s, time domain step is ∆t ≈ 1 min, and the spatial step should be the order of ∆x ≈ 1 mm. Choosing values for the time and step domains depend on the numerical values of the parameters in the equation to be solved, mostly the diffusion coefficient. Besides satisfying Equation (5), the time and space domain steps should be small enough to ensure the accuracy of the solution. There are no explicitly defined mathematical criteria for this procedure and analysis needs to be performed to ensure the stability of the solution. The domains steps should be arbitrarily small but large enough such that the computation time is acceptable. Then domain steps can be lowered, until there is no significant distinction between results. Another way to check the validity of the numerical solutions is to compare it to the analytical solutions where it is possible.
The initial condition inside the origin is density ρ(z = 0; t = 0) = ρ 0 and along the z-axis the CO 2 density equates to zero ρ(z = 0; t = 0) = 0. Starting from the initial conditions, we can iteratively calculate density in the next step up to an arbitrary point in time.
To determine the values of density at the boundary, the boundary conditions need to be defined as ρ(H, t) = ρ 0 and ρ(0, t) = ρ air . Leakage on the surface can be determined from the flux of the boundary with ground j = −D ∂ρ ∂z z=H . The member 2gρ in the brackets of Equation (3) can be roughly estimated and compared with the pressure gradient. It was already mentioned earlier in this paper that pressure gradient is about 10.5 MPa/km, being equal to 10.5 kPa/m. The density of CO 2 decreases from around 650 kg/m 3 to its atmospheric density of 1.87 kg/m 3 . If one applies the arithmetic mean, a value of around ρ = 320 kg/m 3 is obtained, which leads to 2gρ ≈ 6 kPa/m; this member is not small and cannot be neglected. At this point, one can then utilize the effective value of density to obtain an approximation of the component in brackets (Equation (3)), which is constant and equal to ∇P ≡ ∇p − 2gρ eff . Thereupon, the following equation is obtained: By introducing the constants V = − k ϕµ ∇P and K = D ϕ , Equation (6) becomes the standard convection-diffusion equation for incompressible fluids: If this procedure, which is introduced in this manuscript, can be performed and effective density can be introduced, then the problem can be simplified to the case on incompressible fluid. This opportunity is explored in the present paper. The possibility of a reduction in the problem can lead to simplified models, which can be solved in a simplified manner. Another convenience is that there is an analytical solution of the diffusion-convection equation from incompressible fluid: where function f (z, t) is given in Appendix A-Equation (A2). Further steps are to find the solution of convection-diffusion equation for compressible fluids, perform sensitivity analysis on the density as factor ant try to find effective density for which convection-diffusion equation will give the same result as the solution for compressible fluids.

Results and Discussion
Saline aquifers as locations for CO 2 injections have an advantage in storage capacity that oil and gas storages do not. However, CO 2 mobility can lead to supercritical CO 2 spreading within large distances, as well as leading to leakage in locations far from the injection point. During injection, CO 2 forms a plume which spreads in lateral directions and the evolution of the plume can vary depending on the aquifer properties. While CO 2 is trapped beneath a caprock due to residual trapping, in horizontal caprock formations the plume will spread laterally, the sharp edge of which can be determined by Equation (1). The aquifer's porosity and CO 2 viscosity will be the predominant influences on the plume's evolution, and while other parameters can make a considerable impact, their values do not vary drastically from aquifer to aquifer. Figure 2 depicts the evolution of the plume based on the aquifer's porosity. Permeability is taken to be k = 20 mD, viscosity is µ c = 0.0395 mPa · s and µ w = 0.312 mPa · s, thickness of the brine is H = 30 m, and the injection rate is Q = 1200 m 3 /day [90]. The injection point lies in the origin of the coordinate system, and the CO 2 plume evolution is showcased for up to thirty years since the injection period. With porosity as a parameter, the full curve to the right illustrates ϕ = 0.1. Each curve to the left depicts the increase in porosity by 2%, whereas the left end curve illustrates ϕ = 0.18 [91]. This is the expected range of porosity in aquifers.  Similar results have been reported by Lindeberg [81], wherein brine had been modeled twenty-five years after injection. The brine had expanded around 4 km from the source of the reservoir, with a permeability of 2000 mD. After 300 years, the brine had reached 7 km, but after a period of 10,000 years, Lindeberg [81] reported distances larger than 10 km, with a permeability of 250 mD. Therefore, porosity has a large impact on Similar results have been reported by Lindeberg [81], wherein brine had been modeled twenty-five years after injection. The brine had expanded around 4 km from the source of the reservoir, with a permeability of 2000 mD. After 300 years, the brine had reached 7 km, but after a period of 10,000 years, Lindeberg [81] reported distances larger than 10 km, with a permeability of 250 mD. Therefore, porosity has a large impact on plume evolution. A porosity variation of 8% influences the radial expansion of the plume by 30%. It is expected that after 30 years, CO 2 can also migrate up to 10 km from its injection point.
Another parameter which has a considerable impact on CO 2 migration is the viscosity of CO 2 . Figure 3 depicts the plume evolution in accordance with parameters in Figure 2, wherein the porosity is fixed at ϕ = 0.15, but the CO 2 viscosity varies in the range Viscosity clearly has an influence over plume evolution; the plume's shape and migration range depend more on viscosity than on porosity. A plume can reach distances of up to 10 km from the injection point, and since CO2 sequestration is planned long-term, a thirty-year period is relatively short for the time needed to fully trap CO2. In order to ensure that dissolution and mineralization, both of which are lengthy processes, complete the CO2 trapping, CO2 should remain trapped for longer periods of time.
Assuming that CO2 injection inside the aquifer occurs in accordance with the abovementioned parameters and stops after thirty years, the plume will evolve until it reaches an equilibrium state. Since the injection term is equal to 0 ( 0 Q  ), the governing equations will need to be modified. We will presently ignore the progress of its evolution into a steady state and try to assess the maximum radial distance of the CO2 migration. Assuming that CO2 forms a relatively thick bulb beneath the bedrock, the evaluation can be simplified; the thickness will vary with distance but in a steady state the variations will be minimal. After the injection period, the plume will reach an arbitrary state. In case of horizontal saline aquifers with structural trapping, the maximum distance from the injection site will be estimated from the effective thickness of the plume: Here, heff is the effective plume thickness which is equal to the mean thickness of the plume in its equilibrium state, long after the injection period [81]. Radial migration will eventually be stopped due to the decreasing of hydrostatic vertical pressure and bouncy force, which depend on the plume thickness. The diffusion and convection of the plume into water will minimize spreading and CO2 will reach its maximum radial distance from the source, which will be in dynamic equilibrium with all processes occurring in the plume and the interface with brine [83]. Figure 4 shows how the maximum radial distance varies depending on the effective Viscosity clearly has an influence over plume evolution; the plume's shape and migration range depend more on viscosity than on porosity. A plume can reach distances of up to 10 km from the injection point, and since CO 2 sequestration is planned long-term, a thirty-year period is relatively short for the time needed to fully trap CO 2 . In order to ensure that dissolution and mineralization, both of which are lengthy processes, complete the CO 2 trapping, CO 2 should remain trapped for longer periods of time.
Assuming that CO 2 injection inside the aquifer occurs in accordance with the abovementioned parameters and stops after thirty years, the plume will evolve until it reaches an equilibrium state. Since the injection term is equal to 0 (Q = 0), the governing equations will need to be modified. We will presently ignore the progress of its evolution into a steady state and try to assess the maximum radial distance of the CO 2 migration. Assuming that CO 2 forms a relatively thick bulb beneath the bedrock, the evaluation can be simplified; the thickness will vary with distance but in a steady state the variations will be minimal. After the injection period, the plume will reach an arbitrary state. In case of horizontal saline aquifers with structural trapping, the maximum distance from the injection site will be estimated from the effective thickness of the plume: Here, h eff is the effective plume thickness which is equal to the mean thickness of the plume in its equilibrium state, long after the injection period [81]. Radial migration will eventually be stopped due to the decreasing of hydrostatic vertical pressure and bouncy force, which depend on the plume thickness. The diffusion and convection of the plume into water will minimize spreading and CO 2 will reach its maximum radial distance from the source, which will be in dynamic equilibrium with all processes occurring in the plume and the interface with brine [83]. Figure 4 shows how the maximum radial distance varies depending on the effective thickness of the CO 2 bulb. The results show that the plume can reach radial distances up to tens of kilometers from the injection point. This result differs drastically from the estimated evolution of plumes during relatively short time period of CO 2 injection. Thirty years of an injection period is relatively small compared to much larger periods of time for plume evaluation.  , the radial transport speed amounts to around 2.6 m a year, showing that it takes 3.8 millennia to reach the radial distance of 10 km.
If the plume provides an opening for CO2 to migrate upwards, the governing equation of upward migration is Equation (2), whose solution was found by employing the described numerical methods expressed with Equation (4). The velocity of the upward flow depends on the density of CO2 and the transport equation might be reduced to the convection-diffusion equation for incompressible fluids by introducing an effective vertical density.
Let us at this point assume that the vertical migration is driven only by diffusion without the convection process in semi-infinite media by setting the superficial velocity equal to zero. Density in that case can be calculated as: In our problem, upward motion is bounded at the end with the surface of the ground, where CO2 enters the atmosphere. The newer, less evolution of diffusion in semi-infinite media will be the same as in the bounded media, until CO2 reaches surface. Afterwards, the diffusion will not be the same but can enable an estimation of the vertical breach time. From Figure 5, it can be seen that the vertical breaches for saline aquifer at the depth of 700 m can occur after about 50 millennia. Figure 5 is based on the Equation (10) for unit 3 2  Figure 4. Dependence of maximum radial range extent of the plume on effective plume thickness during its equilibrium state thirty years after injection, with volumetric injection rate at Q = 1200 m 3 /day [90].
The durations of CO 2 migration and reaching the maximum range vary significantly in relation to CO 2 permeability and viscosity. The flux can be derived as j = − k µ c ∆ρgS [81], wherein ∆ρ presents the difference in density between two phases and S is the upward slope. If the permeability is k = 2000 mD, µ c = 0.06 mPa · s, and h = 0.1 m, leading to ∆ρ = 250 kg/m 3 and the upward slope is S = 10 −3 , the radial transport speed amounts to around 2.6 m a year, showing that it takes 3.8 millennia to reach the radial distance of 10 km.
If the plume provides an opening for CO 2 to migrate upwards, the governing equation of upward migration is Equation (2), whose solution was found by employing the described numerical methods expressed with Equation (4). The velocity of the upward flow depends on the density of CO 2 and the transport equation might be reduced to the convectiondiffusion equation for incompressible fluids by introducing an effective vertical density.
Let us at this point assume that the vertical migration is driven only by diffusion without the convection process in semi-infinite media by setting the superficial velocity equal to zero. Density in that case can be calculated as: In our problem, upward motion is bounded at the end with the surface of the ground, where CO 2 enters the atmosphere. The newer, less evolution of diffusion in semi-infinite media will be the same as in the bounded media, until CO 2 reaches surface. Afterwards, the diffusion will not be the same but can enable an estimation of the vertical breach time. From Figure 5, it can be seen that the vertical breaches for saline aquifer at the depth of 700 m can occur after about 50 millennia. Figure 5 is based on the Equation (10) for unit initial density and diffusion coefficient D = 0.2 · 10 −3 cm 2 /s.
The time domain step is taken to be 0.1 s, while the spatial domain step is set by Equation (5)    , the effective transport of the incompressible flow is equivalent to the convective transport of compressible media, as depicted in Figure 6. The same can be applied to the flow of various fluids with arbitrary parameters. This is a crucial finding in this work and it enables simplifying complex problems with no loss in generalization. Compressible fluids can be modeled as unoppressive with effective density. This can be quite a breakthrough since there are analytical solutions for the later fluids and their properties and behavior in saline aquifers can be easily studied.
The time domain step is taken to be 0.1 s, while the spatial domain step is set by Equation (5). The sensitivity of the results on the time step value is performed to check the stability of the solution. The time step ∆t should be small enough to provide stability of the solution but taking a small value will have a negative effect as the computation time can be quite high. The authors lowered the value of the time step domain until they reached practically unchangeable results. Further lowering of the time step domain produces relative difference in the results expressed in parts of the percentage. Further in the results, the numerical solutions are compared with the analytical solutions in order to verify the validity of the presented results.
The variation of ρ eff illustrates that the value needs to be higher than average, even the maximum of density ρ 0 , to equate the solution for an incompressible fluid to the solution for compressible ones. By estimating that ρ eff = 1.4ρ 0 , the effective transport of the incompressible flow is equivalent to the convective transport of compressible media, as depicted in Figure 6. The same can be applied to the flow of various fluids with arbitrary parameters. This is a crucial finding in this work and it enables simplifying complex problems with no loss in generalization. Compressible fluids can be modeled as unoppressive with effective density. This can be quite a breakthrough since there are analytical solutions for the later fluids and their properties and behavior in saline aquifers can be easily studied. To emphasize the importance of simplification, it is necessary to use exact (Equation (8)) and numerical (Equation (4) Time is taken as a parameter.
To emphasize the importance of simplification, it is necessary to use exact (Equation (8)) and numerical (Equation (4)) solutions of the convection-diffusion equations. Figure  7 illustrates solutions for parameters The analytical solution is calculated for the first 100 components of the sum in Equations (8) and (A1) of the Appendix A. The equilibrium is achieved 50 s after convectiondiffusion goes into motion, but the analytical solution shows deviations due to oscillations in the solution. Adding more components of the sum would refine the solution but would make computing impractical. While the numerical solution would provide beneficial results for the given timescale, the long evaluation period, which entails millennia, renders this computational method useless. Utilizing the small values of the convective flow and the diffusion coefficient would be impractical and time-consuming. Parallel computing could significantly improve the outcome, but the issue lies with the divergence of numerical solutions when it comes to large scales. The discretization process would also have to be smaller in scale than the diffusion coefficient and the velocity of the fluid flow require it to be. Additionally, a scale does not consider the spatial dimensions of the problem. The condition for numerical solution stability regarding diffusion is given by Equation (5), which gives the constraints on spatial discretization regarding diffusion coefficient. Consequently, simulating periods of millennia are constrained with spatial dimensions. The discretization of space needs to be order of meters in maximum for hundred meters of spatial size. On the other side, the diffusion coefficient constrains time discretization to order of seconds. To simulate millennia second by second is like measuring the distance to the moon with micrometers. Modeling vertical migration in the equilibrium state is gaining attention due to fast and accurate simulation [92].
Taking into consideration prior analyses, the convection flow of the fluid can be described with an average speed of V = 3.5 × 10 −3 m/s and permeability of k = 2000 mD, the value of which, albeit not quite accurate, results in a worst-case scenario, wherein the convective flow brings CO 2 to the surface in the span of two months. Depending on the effective surface area of cracks through which CO 2 can leak, its return into the atmosphere could be slow-moving. Nevertheless, horizontal aquifers, in any given scenario, pose the possibility of CO 2 returning into the atmosphere before mineralization manages to permanently entrap CO 2 .

Conclusions
By mathematical modeling horizontal aquifer, the authors evaluated the storing possibilities of industrial CO 2 . The analysis performed in this work demonstrates that many factors pose potential risks to the process of CO 2 sequestration and its sustainability, with the mobility of CO 2 being one of the main ones; it is of interest to have a higher mobility of CO 2 since it guarantees larger storage capacities but a higher mobility would increase the chances of CO 2 escaping the storage and returning into the atmosphere.
Storage durability depends on the type of aquifer and its characteristics. Knowledge of the evolution of storage processes in a particular aquifer is not deterministic, rather we can talk about estimates and the probability of some outcomes in the sequestration process. It is certain, however, that the porosity of the material forming the aquifer, as well as the CO 2 permeability, have a considerable impact on CO 2 migration. The porosity of aquifer brines can vary from 10% to about 20%, and for a given range the evolution of plume and its radial distance can vary up to 30%. This is an important factor because it determines mobility and capacity at the same time in inverse proportionality. Saline water viscosity is the second most important parameter and can vary in range (0.1 − 0.4) mPa/s. For a given variation of water viscosity, the radial distance of plume can be doubled in range from 5 to 10 km for 30 years injection period. In such a way, the formed brine after injection period will continue to evolve and have diffusion and convection migrations until reaching some equilibrium state. This research showed that, depending on the parameters of the aquifer, brine can evolve into a stadium for which the maximal radial distance can be up to several tens of kilometers from the injection point in worst-case scenarios for injection periods of 30 years. For long injection periods and large injection rates, the radial migration of CO 2 can be even larger.
Long injection periods and high injection rates can, however, increase the radial migration of CO 2 and leakage becomes probable when CO 2 finds an upward path inside the aquifer. The results of this work showed that in the worst-case scenarios, where breaches occurred at open or poorly sealed abundant wells, leakage of the CO 2 can occur after a very short period, which can be measured in months and years. For millennium-lasting projects, this is quite a small time.
Poor seals can have high permeability, and in the long-term trapping of CO 2 it will allow considerable brine leakage out of the formation vertically upwards. The paper has The analytical solution of the diffusion-convection equation for incompressible fluid, given by Equation (8), consists of function f (z, t), which has the form: