Numerical Prediction of the Behavior of CO 2 Bubbles Leaked from Seaﬂoor and Their Convection and Di ﬀ usion near Southeastern Coast of Korea

: Among various carbon capture and storage technologies to mitigate global warming and ocean acidiﬁcation due to greenhouse gases, ocean geological storage is considered the most feasible for Korea due to insu ﬃ cient inland space to store CO 2 . However, the risk of CO 2 leakage and the behavior and environmental e ﬀ ects of the leaked CO 2 need to be assessed for its successful implementation. Therefore, the behavior of CO 2 bubbles / droplets dissolving into the surrounding seawater and the di ﬀ usion of dissolved CO 2 by ocean ﬂows should be accurately predicted. However, ﬁnding corresponding research has been di ﬃ cult in Korea. Herein, the behavior and convection-di ﬀ usion of CO 2 that was assumed to have leaked from the seaﬂoor near the southeastern coast of Korea were numerically predicted using a multi-scale ocean model for the ﬁrst time. In the simulation region, one of the pilot projects of CO 2 ocean geological storage had started but has been temporarily halted. In the ocean model, hydrostatic approximation and the Eulerian–Lagrangian two-phase model were applied for meso- and small-scale regions, respectively. Parameters for the simulations were the leakage rate and the initial diameter of CO 2 . Results revealed that all leaked and rising CO 2 bubbles were dissolved into the seawater before reaching the free surface; further, the change in the partial pressure of CO 2 did not exceed 500 ppm during 30 days of leakage for all cases.


Introduction
Carbon capture and storage (CCS) is one of the technologies used to mitigate global warming and ocean acidification. Various CCS methods are available and can be categorized into inland or ocean geological storage, direct injection and so on. Among these, only inland and ocean geological storage are considered viable following the prohibition of the direct injection method due to its uncertain environmental impacts. CO 2 geological storage is a method for capturing CO 2 from power plants or industrial processes without releasing it into the atmosphere, transporting it to sites suitable for geological storage, and storing it stably underground for the long term. Furthermore, it has the advantage of large-scale CO 2 reduction. Geological storage has been mainly carried out in countries with sufficient inland storage space. However, for countries with limited inland space such as Korea, Japan, Norway, and the United Kingdom, ocean geological storage is considered the most feasible option.
The Sleipner CO 2 injection project in Norway was the world's first industrial offshore CO 2 ocean geological storage project, through which more than 16 Mt of CO 2 was injected from 1996 to 2014. The injection rate was approximately 0.9 Mt/y during the early years and reduced slightly during later years due to reduced gas flow from the Sleipner Vest [1,2]. From 2008 to 2012, 1.6 Mt of CO 2 was injected into the Snøhvit gas field. Although the injection was occasionally halted due to operational challenges at the LNG plant during that period, approximately 0.5 Mt of CO 2 was stored in the well, and injection has continued since 2011 [3]. A large-scale CCS demonstration project in Japan's Tomakomai area, which can store 0.1 Mt/y of CO 2 in two reservoirs that lie 1100 m and 2400 m below the seabed, is being undertaken by the Japanese government [4]. In the case of Korea, a CO 2 storage project that can store 1 Mt/y was being constructed near the southeastern coast of Korea, but was temporarily stopped in 2018.
To utilize any kind of CCS technology, public acceptance must be ensured through the risk assessment of CO 2 leakage, monitoring the behavior of leaked and dissolved CO 2 (DCO 2 ), and environmental assessment of the leaked CO 2 . The Quantifying and Monitoring Potential Ecosystem Impacts of Geological Carbon Storage (QICS) project in the UK was created to assess the environmental impact of CO 2 release experiments in Ardmucknish Bay [5]. Sellami et al. [6] investigated the dynamics of leaked CO 2 bubbles in a plume in the bay through observational data obtained from the QICS project. Because of the limitations of field experiments, studies utilizing computational fluid dynamics (CFD) have been carried out to estimate the behaviors of leaked CO 2 and DCO 2 . The highly complex phenomena that occur when liquid CO 2 is injected into the deep ocean was numerically studied by Eulerian-Eulerian two-phase CFD simulations by Alendal and Drange [7] and Eulerian-Lagrangian simulations by Sato and Sato [8]. Jeong et al. [9] developed a multi-scale ocean model, which builds a bridge between smaller near-field scale and larger regional-scale models to predict the fate of DCO 2 in the deep ocean. Kano et al. [10] conducted numerical simulations on the behavior of CO 2 bubbles and droplets leaked from the seafloor into water columns in uniform flows perpendicular to the leakage band. They showed that CO 2 dissolved in seawater before returning to the air under their simulation conditions, which indicates that the ocean can play the role of a buffer that does not allow CO 2 to return to the atmosphere. Kano et al. [11] developed a multi-scale numerical method to predict the behavior of dissolvable CO 2 bubbles leaked from the seafloor and dissolved mass in the ocean. A simulation using this model was conducted with real topography and tidal currents near the Japanese coastline. Mori et al. [12] used the model developed by Kano et al. [11] and conducted case studies in which the ratio of CO 2 seepage in the dissolved phase and the proportion remaining in the sediment were changed to predict CO 2 concentration distributions in Ardmucknish Bay. However, in Korea, finding corresponding research has been difficult, and the few studies, such as that by Kang et al. [13], have only focused on the behavior inside CO 2 reservoirs and interactions between the reservoir and the surface of seabed through the fault. Such studies cannot be applied to consider ocean flows nor assess environmental impacts in the ocean.
In this study, the behavior and convection-diffusion of CO 2 , which was assumed to have leaked from the seafloor near the southeastern coast of Korea, where a pilot project of CO 2 ocean geological storage started but was temporarily stopped, were numerically predicted using a multi-scale ocean model for the first time. In the ocean model, hydrostatic approximation and Eulerian-Lagrangian two-phase models were applied for the meso-and small-scale regions, respectively. Numerical simulations involving changes in the main parameters, such as the initial diameter of CO 2 and leakage rate, were carried out to investigate its effects on the simulation results, especially regarding environmental impacts.

Methodology
The multi-scale ocean model used in this study is an improved version of the original Maritime Environment Committee (MEC) ocean model developed by the Japan Society of Naval Architecture and Ocean Engineers (JASNAOE). In the MEC model, hydrostatic approximation and non-hydrostatic Appl. Sci. 2020, 10, 4237 3 of 16 (i.e., full-3D) models are applied for the meso-and small-scale domains, respectively. In the mesoscale domain, tidal flow is generated under hydrostatic (pressure) approximation to reduce the computational time. The applicability and accuracy of the hydrostatic model of the MEC model can be found in the research of Kano et al. [11], Lee et al. [14], and so on. A full-3D small-scale model is necessary when the vertical flow component cannot be ignored. The spatial connection between the two models is such that the full-3D model domain matches one grid-column of the hydrostatic model. Kano et al. [11] modified the MEC model by adopting a two-phase model for the small-scale domain, where the continuous liquid phase (seawater) and dispersed gas/liquid phase (individual CO 2 bubble/droplet) are solved using Eulerian and Lagrangian methods, respectively. This model, hereafter termed the "MEC-CO 2 model", was used in the present study.
The governing equations of the mesoscale model adopting the hydrostatic approximation are as follows.
∂u ∂x where u, v, and w are the velocity components in the x, y, and z directions, respectively. p, ρ, g, and f are the pressure, density of seawater, gravity acceleration, and Coriolis force coefficient, respectively. A M and K M are the horizontal and vertical eddy viscosities, respectively. The transport equations of scalar properties φ (i.e., temperature (T), salinity (S), and DCO 2 (C)) are as follows: where A D and K D are the horizontal and vertical eddy diffusivities, respectively. As mentioned earlier, the Eulerian-Lagrangian two-phase model is applied for small-scale full-3D regions. The continuity and Navier-Stokes equations for the continuous phase are as follows: where, where α is the volume fraction and u is the velocity vector. n d , V, and V cell are the number of bubbles in a computational cell, volume, and the volume of a computational cell, respectively. Γ is the mass transfer from a bubble at the interface and P is the pressure. κ and ν are the thermal diffusivity and kinematic viscosity, respectively. The subscripts c and d denote the continuous water phase and dispersed CO 2 bubble phase, respectively. ρ d and ρ c are given by Pitzer and Sterner [15] and Alendal and Drange [7], respectively, and ν t is determined by Smagorinsky's model [16]. The transport equations of T, S, and C are also solved as given below: where C I and C cell are C at the bubble surface and in a computational cell, respectively. D is the diffusion coefficient. Pr t and Sc t are the turbulent Prandtl and Schmidt numbers, respectively. Hirai et al. [17] found that C I matches the solubility of CO 2 , which was measured by Weiss [18]. The mass transfer rate k is modeled as in Chen et al. [19] and d e is the equivalent diameter of a bubble. For the dispersed phase, the mass conservation and motion equation are solved for each bubble in the Lagrangian frame as: where, where β is the coefficient of the added mass of a bubble and g is the gravity vector. f D , f L , and C D , C L are the drag and lift forces and the coefficients of a bubble, respectively. u r is the relative velocity of the dispersed phase to the continuous phase, and ω c is the vorticity of a bubble. β and C L are both set to 0.5. The details on the modeling C D of a bubble can be found in [19]. The governing equations were discretized by the finite volume method using orthogonal and staggered grids: velocity components were defined at the cell faces, and the other variables were defined at the cell centers. For spatial discretization, third-order up-winding and second-order central differencing schemes were adopted for the advection and diffusion terms, respectively. The second-order Adams-Bashforth method was used for the time integration in the full-3D model. Figure 1 shows the target area of the present simulations near the southeastern coast of Korea, where one of the candidate sites for CO 2 ocean geological storage is located. The data for inland and seafloor topographies were obtained from the Global 30 Arc-Second Elevation Dataset (GTOPO30) of the United States Geological Survey's Earth Resources Observation and Science Center and the JODC-Expert Grid Data for Geography-500m (J-EGG500) of the Japan Oceanographic Data Center, respectively. Grid systems for the mesoscale domain were generated using the preprocessor of the MEC model with these data. The dimensions of the computational domain were approximately 100 km × 50 km × 1 km, and the number of grid cells were 50 × 50 × 63 in the x, y, and z directions, respectively, as shown in Figure 2a. The vertical cell column contains the full-3D small-scale domain, of which the size is 2 km × 1 km × 0.18 km with 100 × 100 × 28 grid cells and is shown by a closed cell in Figure 2b. The vertical grid levels are listed in Table 1, where layers 1 and 28 match the free surface and the depth of CO 2 leakage, that is, 190 m, respectively.

Simulation Conditions
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 16 JODC-Expert Grid Data for Geography-500m (J-EGG500) of the Japan Oceanographic Data Center, respectively. Grid systems for the mesoscale domain were generated using the preprocessor of the MEC model with these data. The dimensions of the computational domain were approximately 100 km × 50 km × 1 km, and the number of grid cells were 50 × 50 × 63 in the x, y, and z directions, respectively, as shown in Figure 2a. The vertical cell column contains the full-3D small-scale domain, of which the size is 2 km × 1 km × 0.18 km with 100 × 100 × 28 grid cells and is shown by a closed cell in Figure 2b. The vertical grid levels are listed in Table 1, where layers 1 and 28 match the free surface and the depth of CO2 leakage, that is, 190 m, respectively.  To reproduce proper ocean flows, the major tidal components of M2, O1, K1, and S2, which were obtained from the NAO99b model [20] and listed in Table 2, were interpolated and imposed on the Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 16 JODC-Expert Grid Data for Geography-500m (J-EGG500) of the Japan Oceanographic Data Center, respectively. Grid systems for the mesoscale domain were generated using the preprocessor of the MEC model with these data. The dimensions of the computational domain were approximately 100 km × 50 km × 1 km, and the number of grid cells were 50 × 50 × 63 in the x, y, and z directions, respectively, as shown in Figure 2a. The vertical cell column contains the full-3D small-scale domain, of which the size is 2 km × 1 km × 0.18 km with 100 × 100 × 28 grid cells and is shown by a closed cell in Figure 2b. The vertical grid levels are listed in Table 1, where layers 1 and 28 match the free surface and the depth of CO2 leakage, that is, 190 m, respectively.  To reproduce proper ocean flows, the major tidal components of M2, O1, K1, and S2, which were obtained from the NAO99b model [20] and listed in Table 2, were interpolated and imposed on the  To reproduce proper ocean flows, the major tidal components of M2, O1, K1, and S2, which were obtained from the NAO99b model [20] and listed in Table 2, were interpolated and imposed on the open boundaries of the mesoscale model domain with nonreflecting boundary conditions by Hino and Nakaze [21].  open boundaries of the mesoscale model domain with nonreflecting boundary conditions by Hino and Nakaze [21].  Figure 3 shows the initial conditions of T, S, and DCO2 obtained from the data of the Array for Real-time Geostrophic Oceanography project by the National Institute of Meteorological Sciences of Korea (http://argo.nims.go.kr). To calculate the fluxes of temperature and salinity at the free surface, climate data, which were the temporal average of observed data by the National Aeronautics and Space Administration (NASA, https://power.larc.nasa.gov/data-access-viewer/), were used, as depicted in Table 3.  To calculate the fluxes of temperature and salinity at the free surface, climate data, which were the temporal average of observed data by the National Aeronautics and Space Administration (NASA, https://power.larc.nasa.gov/data-access-viewer/), were used, as depicted in Table 3. The leakage rate and initial diameter of the CO 2 bubble were selected as the main parameters. The assumed leakage rates for the present study were 3800, 50,000, and 100,000 t/y, based on the study of Kano et al. [11], where two cases for the leakage rate were studied: an extreme case, 94,600 t/y, which assumed that a large fault accidentally connected the CO 2 reservoir and the seafloor; and a reasonable case, 3800 t/y, based on the seepage rate of an existing enhanced oil recovery (EOR) site [22]. The initial diameters of the CO 2 bubble were set to 5, 10, and 20 mm. The simulation cases and leakage are listed and illustrated in Table 4 and Figure 4, respectively.  [22]. The initial diameters of the CO2 bubble were set to 5, 10, and 20 mm. The simulation cases and leakage are listed and illustrated in Table 4 and Figure 4, respectively.  Quantitative criteria are needed to assess the environmental impact when CO2 leakage occurs. Kikkawa et al. [23] elucidated that the biological impacts of CO2 in the ocean should not be related solely to pH, but also to the partial pressure of CO2 (pCO2). Kita and Watanabe [24] collected LC50 and LT50 data against pCO2 for various marine species and proposed that the change in pCO2 (ΔpCO2) of 5000 ppm is the no-observed effect concentration (NOEC) and that 500 ppm is the predicted no-effect concentration (PNEC). Although these values have not been authorized, we refer to them as tentative standards in this study.  Quantitative criteria are needed to assess the environmental impact when CO 2 leakage occurs. Kikkawa et al. [23] elucidated that the biological impacts of CO 2 in the ocean should not be related solely to pH, but also to the partial pressure of CO 2 (pCO 2 ). Kita and Watanabe [24] collected LC50 and LT50 data against pCO 2 for various marine species and proposed that the change in pCO 2 (∆pCO 2 ) of 5000 ppm is the no-observed effect concentration (NOEC) and that 500 ppm is the predicted no-effect concentration (PNEC). Although these values have not been authorized, we refer to them as tentative standards in this study.   Figure 7 shows contour maps of some variables on the center plane in the x direction (x = 1000 m) of the small-scale domain, 30 days after the start of leakage. From the distribution of void rate and cell-averaged number densities of CO2 bubbles, it can be seen that all CO2 bubbles leaking from the depth of 190 m would be dissolved at a depth of 70 m before reaching the free surface. Because the initial diameter of the CO2 bubbles are the same, but the leakage rates are different for these cases, the void rates of the bubbles near the leakage depth are higher than other depths and become lower as the depth becomes shallower due to the rising and dissolution of the bubbles to the surrounding seawater. The reason for the higher values of the number densities near the depth of 100 m is that the density of the bubbles is almost the same as that of seawater, and the volume of the bubbles is very small, which makes the buoyant force zero. Although the velocity vector fields or contour maps of the velocities are not presented, one may know that complex flows exist owing to the rising of the CO2 bubbles with surrounding seawater and tidal flows from the distribution of the DCO2 shown in Figure 7c. Since more bubbles exist and dissolution lasts for a longer time as the leakage rate increases, relatively higher values and wider areas of dispersed DCO2 are observed, as shown in the figure. Comparing the contour maps of DCO2 with those of the changes in DCO2 (ΔDCO2) due to the dissolution of the bubbles, no large discrepancy is observed because the background DCO2 is much smaller.  Figure 7 shows contour maps of some variables on the center plane in the x direction (x = 1000 m) of the small-scale domain, 30 days after the start of leakage. From the distribution of void rate and cell-averaged number densities of CO 2 bubbles, it can be seen that all CO 2 bubbles leaking from the depth of 190 m would be dissolved at a depth of 70 m before reaching the free surface. Because the initial diameter of the CO 2 bubbles are the same, but the leakage rates are different for these cases, the void rates of the bubbles near the leakage depth are higher than other depths and become lower as the depth becomes shallower due to the rising and dissolution of the bubbles to the surrounding seawater. The reason for the higher values of the number densities near the depth of 100 m is that the density of the bubbles is almost the same as that of seawater, and the volume of the bubbles is very small, which makes the buoyant force zero. Although the velocity vector fields or contour maps of the velocities are not presented, one may know that complex flows exist owing to the rising of the CO 2 bubbles with surrounding seawater and tidal flows from the distribution of the DCO 2 shown in Figure 7c. Since more bubbles exist and dissolution lasts for a longer time as the leakage rate increases, relatively higher values and wider areas of dispersed DCO 2 are observed, as shown in the figure. Comparing the contour maps of DCO 2 with those of the changes in DCO 2 (∆DCO 2 ) due to the dissolution of the bubbles, no large discrepancy is observed because the background DCO 2 is much smaller.

Effect of Leakage Rate (Cases 1, 2, and 3/Bubble Size = 20 mm, Leakage Area = 20,000 m)
The contour maps of ∆pCO 2 on the center planes of the x, y, and z directions of the small-scale domain are depicted in Figure 8, where not only the diffusion but also the advection of ∆pCO 2 by the tide are clearly shown. A high ∆pCO 2 is observed near the leakage depth, and the maximum ∆pCO 2 is lower than 200 ppm and does not exceed the PNEC. It can also be seen that the distributions of ∆pCO 2 shown in Figure 8a are similar to those of ∆DCO 2 on the same plane (Figure 7c).
The contour maps of ∆pCO 2 on the constant z planes of the mesoscale domain with respect to depths and elapsed time after leakage for the most extreme case are illustrated in Figures 9 and 10, respectively, where the black rectangular box indicates the full-3D small-scale domain including the leakage area. Both figures show that high ∆pCO 2 regions are placed at depths between 150 and 170 m, not in the leakage depth. This might be due to the convection of DCO 2 by the geostrophic flows. It can also be found that diffusion is more dominant than convection as water depth becomes deeper. The time change of maximum ∆pCO 2 in the mesoscale domain is not significant, as shown in Figure 10. The maximum ∆pCO 2 is lower than 50 ppm and does not exceed the PNEC. Leakage rate = 3800t/y Leakage rate = 50,000t/y Leakage rate = 100,000t/y    Contour maps of some variables on the center plane in the x direction (x = 1000 m) of the small-scale domain at 30 days after the start of leakage are shown in Figure 11. As seen in the distribution of void rate and cell-averaged number densities of CO 2 bubbles, as the initial diameter of the bubbles increases, the rising distance of the bubbles increases, which results from the larger buoyant force of the relatively large volume of the bubble. Because the leakage rates are the same for these cases, the number of bubbles is high within the short band when the initial diameter is 5 mm. Comparing the distributions of DCO 2 or ∆DCO 2 among the cases, one can find that the distribution when the diameter is 10 mm is quite different from those of 5 and 20 mm. The reason for this phenomenon is that the mass transfer rate of a CO 2 bubble is largely affected by the shape of the bubble. The total mass transfer per unit volume decreases as the diameter of a bubble decreases. However, if the diameter exceeds a certain value and changes to the spherical cap, there is a sudden jump, and it decreases as the diameter increases. The threshold is approximately 18 mm. Therefore, the mass transfer rates when the diameters are 5 and 20 mm are similar, but the rate becomes smaller when the diameter is 10 mm, as discussed by Kano et al. (2009). The lower dissolution of the bubbles results in a low DCO 2 Figure 12 shows the contour maps of ∆pCO 2 on the center planes of the x, y, and z directions of the small-scale domain. The maximum ∆pCO 2 for cases 5 and 3 are similar but do not exceed 500 ppm, which is also clearly seen in Figure 13, where those on the constant x, y, and z planes of the mesoscale domain are illustrated.

Conclusions
In this study, the behavior and diffusion of CO2, assumed to have leaked from the seafloor, were numerically predicted for the first time near the one of the candidate sites for CO2 ocean geological storage in Korea. The behavior of leaked CO2 bubbles were analyzed numerically using a multi-scale ocean model to assess its environmental impacts, which is essential for obtaining public acceptance. The main parameters chosen for the simulation were leakage rates and initial diameters of CO2 bubbles. The former was assumed to be 3860, 50,000, and 100,000 t/y and the latter to be 5, 10, and 20 mm, respectively. A total of five simulations were carried out by combining these parameters.
From the simulation results, it was found that all CO2 bubbles were dissolved into the seawater before reaching the free surface in all cases. As the leakage rate increased, relatively higher values were concentrated near the leakage depth, and wider areas of ΔpCO2 were observed because more bubbles existed and dissolution lasted for a longer time. As the initial diameter of the bubbles increased, the rising distance of the bubbles also increased due to the larger buoyant force of the relatively large volume of the bubble. However, small ΔpCO2 values were estimated when the diameter was 10 mm. For the mass transfer rates, though they were similar for the 5 and 20 mm diameters, they were reduced when the diameter was 10 mm.

Conclusions
In this study, the behavior and diffusion of CO 2 , assumed to have leaked from the seafloor, were numerically predicted for the first time near the one of the candidate sites for CO 2 ocean geological storage in Korea. The behavior of leaked CO 2 bubbles were analyzed numerically using a multi-scale ocean model to assess its environmental impacts, which is essential for obtaining public acceptance. The main parameters chosen for the simulation were leakage rates and initial diameters of CO 2 bubbles. The former was assumed to be 3800, 50,000, and 100,000 t/y and the latter to be 5, 10, and 20 mm, respectively. A total of five simulations were carried out by combining these parameters.
From the simulation results, it was found that all CO 2 bubbles were dissolved into the seawater before reaching the free surface in all cases. As the leakage rate increased, relatively higher values were concentrated near the leakage depth, and wider areas of ∆pCO 2 were observed because more bubbles existed and dissolution lasted for a longer time. As the initial diameter of the bubbles increased, the rising distance of the bubbles also increased due to the larger buoyant force of the relatively large volume of the bubble. However, small ∆pCO 2 values were estimated when the diameter was 10 mm. For the mass transfer rates, though they were similar for the 5 and 20 mm diameters, they were reduced when the diameter was 10 mm.
The estimated maximum ∆pCO 2 for the extreme case was approximately 200 ppm during 30 days of leakage under the present simulation conditions. Therefore, it can be said that the environmental impact caused by the leakage of CO 2 was not significant. The present numerical model is expected to be useful and applicable not only for estimation of the behaviors and environmental impact of leaked CO 2 but for various purposes, such as determination of the proper location for leakage monitoring devices during ocean geological storage of CO 2 .
Scenario-based simulations and parametric studies will be performed in the near future. The main parameters, such as the leakage depth closely related to the phase of the leaked CO 2 , the vertical distributions of the scalar properties according to the seasonal change, and the shape of the leakage area, will be considered. Furthermore, validation and improvement of adopted sub-models will be continued as uncertainties and assumptions remain, although well-known or validated ones were implemented in the present model.