International Journal of Greenhouse Gas Control 110 (2021) 1034412ecosystem (which has been shown to modulate the impacts of climate change through CO2 uptake) that needs to be maintained (Artioli et al., 2014, Barange et al., 2011). The practical adoption of CCS depends on the ability to demonstrate storage integrity in the context of both regulation (Dixon et al., 2015) and public perception (Bradbury, 2019), with the need for a trusted monitoring strategy. A large concern in the past has been storage integrity and the likelihood of leakage (Inter-governmental Panel on Climate Change 2019, Monastersky, 2013), with the corresponding environmental impacts and the need for both storage monitoring strategies and leakage quantification (Dixon et al., 2015). Whilst now accepting that properly engineered storage is not expected to leak, regulatory agencies require assessments of the risk, with the potential impact to be both understood and quantified should CO2 be released into the water column (Blackford et al., 2014, Blackford et al., 2020). Any acute impacts on marine biology would occur within the near field (Shirayama and Thornton, 2005), which covers a scale of the seawater from meters up to kilometres. Full scale experiments, with large releases are potentially economic and environmentally expensive. Therefore small scale field and laboratory experiments are required (Roberts and Stalker, 2020), such as the Quantifying and Monitoring Potential Ecosystem Impacts of Geological Carbon Storage (QICS) project (Blackford et al., 2014), Measurement, Monitoring and Verifi-cation of CO2 Storage (MMV) and the Strategies for Environmental Monitoring of Marine Carbon Capture and Storage (STEMM-CCS) proj-ect experiments (Dean et al., 2020), with the development of numerical models to further understand the mechanisms of leakage from the seabed into the turbulent seawater, which can be relatively cheaply extrapolated to larger scale situations. This provides extensive data, filling the gaps and uncertainties left from small scale experiments alone whilst minimising economic and environmental cost (Blackford et al., 2020, Blackford et al., 2013, Blackford et al., 2008, Blackford and Gilbert, 2007, Greenwood et al., 2015, Hvidevold et al., 2015, Phelps et al., 2015, Dewar, 2016, Dewar et al., 2015, Dewar et al., 2013, Dewar et al., 2013, Sellami et al., 2015, Cazenave, 2021). In low depth releases (<~500 m) highly buoyant CO2 gas bubbles rise as a plume into the water column, dissolving rapidly as the natural waters are undersaturated in CO2 (Dewar et al., 2013, Chen et al., 2009). This dissolved solution, results in an increase of Dissolved Inorganic Carbon (DIC) (Zeebe and Wolf-Gladrow, 2001), and partial pressure of CO2 (pCO2). The solution also dissociates into carbonate, bicarbonate and hydrogen ions, with the increases in hydrogen ions reducing the pH of the seawater. CO2 solution is slightly denser than the background (Song et al., 2002, Song et al., 2013), sinking towards the seabed (Dewar et al., 2015, Dewar et al., 2013, Dewar et al., 2013), and sub-sequently dispersed, driven by local hydrodynamics, surface winds, tidal mixing, turbulent diffusion, and residual circulation. In terms of the spatial extent of a release, depending on the release rate, it has been shown to only impact a very limited region of up to the order of a hundred meters immediately adjacent to the source, (Blackford et al., 2020, Blackford et al., 2008, Dewar et al., 2013, Dewar et al., 2013, Chen et al., 2005). Existing models have several shortcomings beyond the lack of experimental data required to develop precise modelling systems. There is a shortfall in the number of multi-phase numerical models (models that can simulate more than one fluid) designed to predict leakages from low depth, dissolving, bubbly flow for CO2 leakage scenarios (Dewar et al., 2015, Dewar et al., 2013, Dewar et al., 2013). Further to this, investigating leakage timeframes from days to weeks and beyond re-quires analysis in multiple scales (Blackford et al., 2020), investigating the local impacts (Blackford et al., 2020, Dewar et al., 2015, Dewar et al., 2013, Dewar et al., 2013) in the order of meters and the regional impacts (Mori et al., 2015) based on changes in currents and tide in the order of tens of kilometres. Except for the modelling system by Mori, Kano, Sato et al. (Mori et al., 2015, Kano et al., 2010), there are a lack of multi-phase, multi-scale modelling systems allowing data transfer between each scale. Even so, the rectilinear mesh used therein can be computationally inefficient if applied over larger coastal areas. This suggests the need for development of a multi-scale, unstructured mesh model (Chen et al., 2003) for good simulations of turbulent flows around the coastline capable of calculating multi-phase hydrodynamics across the mesh. In this study, a 3D physicochemical modelling system is developed for CO2 leakage analysis in coastal waters, built on the unstructured-grid Finite-Volume Community Ocean Model (FVCOM) (Chen et al., 2003) to provide the local and regional hydrodynamics. The Predicting Leakage Using Multi-phase Equations (PLUME) model is developed and inte-grated as a two-way plume modelling system allowing investigation of localised physicochemical impacts utilising the carbonate system from European Regional Seas Ecosystem Model (ERSEM) (Artioli et al., 2012) to analyse the ecosystem impacts. The modelling system and domain are set on the tidally dominated north-west European continental shelf, utilising the Scottish shelf model mesh and boundary forcing system (CH2M - 2016), providing inputs to internally nested domains for the specific sites for high resolution (500 m to 50 cm) analysis. The high-resolution domain allows localised physical and chemical analysis for both the coastal QICS CO2 release experiment (Blackford et al., 2014) in a semi-enclosed bay off the West coast of Scotland, and the off-shore STEMM-CCS CO2 release experiment in the vicinity of the Goldeneye potential storage complex in the Northern North Sea. From the QICS experiment, the first of its kind, large anomalies were found to occur between the experimental readings and the limited nu-merical models (Dewar et al., 2015). With the assistance of the new numerical modelling system, this study aims to resolve these anomalies in detectability and environmental impact findings and validate findings from both controlled sub-seabed release experiments of CO2 (QICS and STEMM-CCS). In-situ Experiments This section outlines the background to this study, with experimental data and findings. The two experiments have similarities and significant differences as described in Table 1. As the QICS experiment was located in a semi-enclosed bay, the horizontal tidal currents were quite strong, up to 40.5 cm/s (Atamanchuk et al., 2015), and a water depth of ~10 m with a fully mixed water column. This compares with residual currents of up to 25 cm/s in the STEMM-CCS experiment (Flohr et al., 2021) near to the Goldeneye complex at ~100 m depth, giving contrasting hydro-dynamics and a stratified water column. Both releases observe bubbly flow throughout the experiment, where the bubble formation is controlled by the structure of the seabed sedi-ments and the fluid interactions as described in Dewar et al. (Dewar et al., 2015). Higher rates with flow from the seabed under pressure, would require jet dynamics affecting bubble sizes and water hydrodynamics. The QICS Experiment QICS was a project with the aim of improving the understanding of potential impacts that a leakage from carbon dioxide geological storage could have on the ecosystem, and to investigate a variety of techniques and methods that may be suitable for monitoring leakage. The Table 1 The Experiment Characteristics QICS STEMM-CCS Water depth (m) 9-12 ~ 120 Injection depth below the seafloor (m) 15 3 Injection rates (kg/day) 0 - 208 0 – 143 Sediment type sandy Clay M. Dewar et al.
International Journal of Greenhouse Gas Control 110 (2021) 1034413experiment involved drilling a narrow borehole from land, terminating in unconsolidated sediments 10 m below the sea floor, with 9–12 m head of seawater in a semi-enclosed bay as shown in Figure 1a. The release was carried out in May – June 2012, continuing for 37 days with a cu-mulative release of 4.2 tonnes of CO2 (Blackford et al., 2014). The data required to model the release experiment, other than the site-specific data (location, depth, salinity, temperature, currents), are the release rates to the seawater, bubble sizes, and pockmark distribution map. The bubble size distribution and pockmark positions are provided by Sellami et al. (Sellami et al., 2015) as shown in Figure 2a and Figure 2b respectively. The size distribution was recorded from a small sample of the experiment as discussed by Sellami et al. (Sellami et al., 2015). However, due to the low flow rates from the seabed 15 meters above the injector, the gas lost its initial injection pressure and only rose through buoyancy. It is therefore possible to predict the size from the mo-mentum, buoyancy and surface tension as described in Dewar et al (Dewar et al., 2015) without any impact from the flow rate. Changes in pockmarks would however cause different sized bubbles to appear. The release rate to the water column has been taken from Berg`es et al. (Berg`es et al., 2015) which showed the release rate in the final stages of experiment. The pattern that varies with the tide (Blackford et al., 2014) has been extrapolated to the start of the experiment using the same in-jection rate / bubble release rate ratio as shown in Figure 2c. Back-ground values for Total Alkalinity (TA) and DIC were measured as 2307 μmol/kg and 2128 μmol/kg respectively (Anon 2021). Anomalies were found to occur between experimental readings (Blackford et al., 2014, Atamanchuk et al., 2015, Kita et al., 2015, Shi-tashima et al., 2015, Blackford et al., 2015, Lichtschlag et al., 2015) and the numerical models (Dewar et al., 2015, Mori et al., 2015, Maeda et al., 2015). The models under predicted physicochemical changes when simulating the physically measured seepage rate, based only on the observed gas flow rate (~15% of the injection rate (Blackford et al., 2014)), providing maximum levels of pCO2 of 443 – 530 μatm (Dewar et al., 2015, Mori et al., 2015). This is compared to the experimental measured values of an average of 390 - 500 μatm, 30 cm above the seabed, but at times increasing up to 1200 - 1250 μatm, with occasional peak values of 1500 - 1600 μatm from observations 3 cm above the seabed (Blackford et al., 2014, Atamanchuk et al., 2015, Kita et al., 2015, Shitashima et al., 2015), showing the large differences over a very small distance. These findings therefore led to some conclusions that the seepage rates for CO2 (gas bubbles and dissolved solution) from the seabed are still largely unknown, with uncertainties of potential CO2 dissolving in the sediments prior to seepage (Dewar et al., 2015, Mori et al., 2015, Maeda et al., 2015). This provides a build-up of DIC being released with the bubbles when the seepage rate increases at low tide, and small bubbles that couldn’t be measured dissolving quickly in the water column (Dewar et al., 2015, Mori et al., 2015). In this paper another mechanism is investigated; the seepage rate measurements are reasonable, but higher model resolution is required to be able to predict the high peaks, with the tidal cycle also having a large effect on the plume dynamics. This mechanism has been discussed (Dean et al., 2020, Maeda et al., 2015) but has yet to be demonstrated in practice. The STEMM-CCS Experiment The STEMM-CCS experiment, a controlled CO2 release in the central North Sea near the Goldeneye platform, dominated by north-south tidal currents, aimed to expand upon knowledge gained from the QICS experiment amongst others (Flohr et al., 2021) as shown in Figure 1b. This experiment was designed to imitate an unintended release of CO2 from a geological CO2 storage site to the seabed. The main objective was to advance the fundamental understanding of the marine environment impacts above a CO2 storage site, also through detecting, characterising and quantifying gaseous and dissolved CO2, provide cost-effective environmental monitoring and leakage quantification techniques (Flohr et al., 2021). The experiment involved inserting a pipe into un-consolidated sediments 3 m below the seabed, with ~120 m head of seawater in open waters. The release was carried out in May 2019, continuing for just over 12 days. As with the QICS experiment, the release rates, bubble sizes, and pockmark distribution map are required for modelling. These parame-ters were recorded through optical, acoustic and gas collection tech-niques through use of remotely operated vehicles (ROV), autonomous underwater vehicles (AUV), optical landers and hydrophone arrays (Flohr et al., 2021). Bubble sizes were predicted by passive acoustics and gas bubble imaging, giving the bubble size distribution shown in Figure 3a (Li et al. 2021), and a map of the bubble streams was collated as shown in Figure 3b. The gas flow rate into the sediments was initially set to 6 kg/d, and then sequentially increased over the duration of the experi-ment to a maximum of 143 kg/d with a cumulative release of 675 kg of CO2 (Flohr et al., 2021) as shown in Figure 3c. As with the QICs experiment, the size distribution is recorded from a small sample of the experiment as discussed by Li et al. (2021). However, due to the low flow rates from the seabed, the flow rate does not impact the bubble size (Dewar et al., 2015). Changes in pockmarks would however cause different size bubbles to appear. Estimates of the release rate to the seawater derived using the gas bubble sampler’s inverted funnel (seep flow rate measurements) showed that the total seepage rates increased from a minimum of ~1.3 kg CO2 d1 (22%) at the lowest injection rate to a maximum of ~70 kg CO2 d1 (48%) at the highest injection rate (Flohr et al., 2021), also shown in Figure 3c, with the rest retained or ejected as a dissolved solution. This is to be compared with the data from eddy covariance techniques that predicted a best estimate of up to 66 ±19% (Koopmans et al., 2021). Figure 1.Graphical representations of the in-situ experimental setup (not to scale). a) The QICS experiment; b) The STEMM-CCS experiment. M. Dewar et al.
International Journal of Greenhouse Gas Control 110 (2021) 1034414Similarly, measurements of benthic gradients of pH using lab-on-chip sensors yielded estimates of 42±17% (Schaap et al., 2021). No bubbles were found to be visible higher than 8m above the seabed (Flohr et al., 2021) showing the fast dissolution rate. pH changes were recorded using lab-on-chip pH sensors, giving a reduction of up to 0.6 units as a peak at ~2.6 m downstream of the plume, with the typical reduction of ~0.3 units during the largest injection rate (143 kg/d) when the tide pushed the flow past the sensor (Schaap et al., 2021). Background values for TA were measured as 2310 μmol/kg, with the corresponding background DIC derived as 2140 μmol/kg (Schaap et al., 2021). Modelling System This section outlines the principal aspects of the study, with the development and application of the hydrodynamic multi-scale and multi-phase modelling system. To investigate the fate and impacts of leakage into the local environment, model coupling and nesting are required. As shown in the experiments in Section 2, the CO2 plume of bubbles manifests locally from centimetres to tens of meters. In order to investigate and predict the plume dynamics, including rise height, horizontal movements, rate of dissolution, and concentration distribu-tion of the CO2 solution, the Predicting Leakage Using Multi-phase Equations (PLUME) model is developed, with the ability to nest into various oceanic hydrodynamic and biogeochemical modelling systems. Interactions of bubbles with seawater are affected greatly by small scale ocean dynamics, such as tides, currents, and turbulence, which, however, develop at scales from global, coastal, regional, down to the local environment. These cascading hydrodynamics are modelled by FVCOM (Chen et al., 2003), and biogeochemical changes to the water column, delineate the overall fate and impacts in the water column. These are obtained from offline linking of the plume and ocean hydro-dynamic models to the carbonate system (Blackford and Gilbert, 2007, Blackford et al., 2004) from ERSEM (Butensch ̈on et al., 2016). An interaction chart is provided in Figure 4 showing how the models (PLUME, FVCOM, and ERSEM) interact with each other within the modelling system, communicating through the physiochemical and biogeochemical variables. PLUME The PLUME model is developed (based on the governing equations from the individual bubble model by Chen et al. (Chen et al., 2009)), that not only resolves the bubble dynamics, but further provides plume dynamics, and provides source terms for changes in mass, momentum and dissolved solutions in the water column to the ocean hydrodynamic and biogeochemical models. The model is a Lagrangian based system, with a set of non-linear ordinary differential equations, which are numerically solved by using the fourth order Runge-Kutta method. PLUME is a general model designed to simulate the interaction dy-namics of gas bubbles, liquid droplets, solid particles, or a combination, with the surrounding fluids dependant on the field of study. In this model, the dynamics of parcels, or collections containing a set number of similar ‘objects’ whether they are bubbles, droplets, or particles (BDPs), are simulated. The number of BDPs in this parcel can be defined and set by the release rate and size distributions as the input. There are no Figure 2.QICS experimental data. a) The bubble size distribution in terms of diameter, top left; b) The bubble stream pockmark locations, top right; c) The leakage rate measured acoustically (solid line) and extrapolated based on injection rate ratio of the measured data (dashed line), with the blue line showing the experiment injection rate and the red point shows the physical leakage rate measurement, bottom. M. Dewar et al.
International Journal of Greenhouse Gas Control 110 (2021) 1034415Figure 3.STEMM-CCS experimental data. a) The bubble size distribution in terms of diameter from 0-20 mm, top left; b) The bubble stream pockmark locations, top right; c) The leakage rate used in the model 50% of the injection rate based on the physical leakage rate measurements (dashed line), with the blue line showing the experiment injection rate and the red points show the physical leakage rate measurements, bottom. Figure 4.The modelling system interaction chart, showing how the models (PLUME, FVCOM, and ERSEM) interact with each other, and which variables are exchanged. M. Dewar et al.
International Journal of Greenhouse Gas Control 110 (2021) 1034416limitations on the number of the BDPs in the parcel, but they require to have similarity in all physical, chemical and biological properties as defined with the given distribution. For leaked or released CO2 in seawater, PLUME interlinks with the hydrodynamic model to bring in the local fluid properties, including the seawater density, currents, temperature, salinity, along with back-ground concentrations of dissolved solutions (such as CO2). PLUME then outputs impacts from the BDPs on momentum and mass from the dis-solved solution back to the hydrodynamic model as shown in Figure 4. Governing equations Consider a collection of objects or BDPs as a parcel, the BDPs in the parcel interact with the surrounding fluids in mass and momentum ex-changes under a thermal equilibrium condition. The total exchange rates of mass and momentum of the BDPs with respect to time can be pre-sented by conservations of mass and momentum (Newton’s 2nd law) as, dmpdt=No ̇mo+mo ̇No(1) dmpur,p)dt=∑nj=1Fj(2) where No and moare the number of objects, and mass (kg) of each BDPs in the parcel, with mp as the total mass of the parcel, ur is the relative velocity vector of BDPs to that of the surrounding fluids (m/s) and Fj are forces acting on the parcel (N). The subscript p indicates the parcel, o the objects or BDPs, and j the forces acting on the BDP. No is initially set for each parcel through the conservation of mass (Eq. 1), assuming a constant number of objects in the parcel (removing the second term). With a given mass release rate, timestep, density and BDP diameter/size/volume provided by a size distribution, the initial number of BDPs in each parcel can be predicted. Two mechanisms can change the total number of BDPs and the distribution in the parcel after release, the breaking up and collision. Details can be found from the previous studies (Dewar et al., 2015, Yao and Morel, 2004). From the field measurements of QICS (Sellami et al., 2015) and STEMM-CCS (Maeda et al., 2015), it is confirmed that both the breaking up and coalescence were rarely observed from these CO2 bubble plumes, and are shown to have a minor effect (Dewar et al., 2015). Therefore, the changes in total bubble number in the parcel are neglected in this study. The mass exchange rate of the BDPs in the parcel can be calculated by the dissolution rate, ̇mo= koAo(CsCb)(3) where ko(m/s) is the effective mass transfer coefficient, Ao (m2) the total surface area of the BDP, Cs the solubility (kg/m3) and Cb the background concentration in surrounding waters (kg/m3). Combining Eq. (1) with Eq. (3) and resolving the geometry (considering each BDP as a sphere to define the equivalent diameter de), the dissolution rate of each BDP in each parcel can be predicted through the equivalent BDP diameter shrinking, at a rate of ̇deo= 1ρo⎛⎝deo3 ̇ρo+2ko(CsCb)⎞⎠(4) where ρo is density (kg/m3) of the BDPs. The first bracketed term on the right-hand side of Eq. (4) gives the effect of compressibility of the BDP to changes in bubble size. The momentum exchange of BDPs in the parcel with surrounding fluids can be identified by the forces acting on each individual BDP from both drag and buoyancy. Adding these forces into equation (2) and resolving for a spherical object, gives the changes in relative velocity of the parcel ur,p (or the individual BDP with an equivalent size of de) as ̇ur,p=γρ(1γ1ρ)g3u2r,p4decd)ur,p ̇momo(5) where, cd is the drag coefficient of each BDP in the surrounding fluids, g is gravity, γρis the density ratio of surrounding fluids to the BDP. Communication with hydrodynamic and ecosystem models The PLUME model communicates through interacting with oceanic hydrodynamic and ecosystem models. As shown in Figure 4. PLUME extracts the 3D currents, pressure, temperature, salinity and density of seawater, along with background concentrations of dissolved solution from the oceanic hydrodynamic model, meanwhile providing the dis-solved mass, the position of the BDP, and the drag force acting on the surrounding seawater as a momentum source. These exchange rates of mass and momentum to the seawater in the oceanic hydrodynamic model are predicted by PLUME by the following source terms, SMass= ̇mPΔtPVHΔtH(6) SMom= ̇ur,p(VPVHΔtPΔtH)(7) where VP and VH are the volumes of the parcel and the hydrodynamic model’s grid cell corresponds with the parcel position, ΔtP and ΔtH are the time steps of the PLUME model and the hydrodynamic model, respectively. ERSEM ERSEM is a comprehensive biogeochemical and ecological model (Butensch ̈on et al., 2016) from which the carbonate system is extracted to predict changes in pH, pCO2 etc. The carbonate system was first introduced to ERSEM by Blackford and Burkill (Blackford and Burkill, 2002) using the HALTAFALL speciation code from Ingri et al. (Ingri et al., 1967). The changes in pH and pCO2 are calculated from the modelled aqueous CO2, salinity and temperature values, along with background values of DIC, TA and total boron (Blackford and Gilbert, 2007, Butensch ̈on et al., 2016). The DIC can be predicted from the addition of the aqueous CO2 to the background DIC measurements. Further details on the ERSEM carbonate system module can be found in Artioli et al. (Artioli et al., 2012) and Butensch ̈on et al. (Butensch ̈on et al., 2016). FVCOM FVCOM is the base model in this system with which the other models interact. It is a unstructured-grid, coastal hydrodynamic circulation model (Chen et al., 2003), with the atmospheric weather forcing through the free surface. The primitive 3D equations of continuity, momentum, energy and tracers (salinity and CO2 solution) are solved to simulate the multiscale dynamics of the hydrodynamic turbulent flows, which are coupled with the PLUME model (Section 3.1) and outputting to ERSEM (Section 3.2). FVCOM has an increasing user base, utilised in a range of various research areas, and has been widely used for simulating hydrodynamics in coastal and shelf waters due to the bathymetry following mesh system and nesting capabilities for high resolution analysis. As a brief and very limited summary, modelling topics have included the stratification and mixing of temperature and salinity (Chen et al., 2007, Ge et al., 2013, Huang, May 2011, Yang and Khangaonkar, 2008, Zheng and Weisberg, May 2012), marine renewables (Cazenave et al., 2016, De Dominicis et al., 2017, Li et al., May 2020), ecosystem and marine habitats (Torres and Uncles, 2011, Aleynik et al., 2016, Waggitt et al., 2018, Waggitt et al., 2016, Waggitt et al., 2016, Adams et al., 2014, Adams et al., M. Dewar et al.