Detection and quantification of a release of carbon dioxide gas at the seafloor using pH eddy covariance and measurements of plume advection

We detected a controlled release of CO 2 (g) with pH eddy covariance. We quantified CO 2 emission using measurements of water velocity and pH in the plume of aqueous CO 2 generated by the bubble streams, and using model predictions of vertical CO 2 dissolution and its dispersion downstream. CO 2 (g) was injected 3 m below the floor of the North Sea at rates of 5.7 – 143 kg d (cid:0) 1 . Instruments were 2.6 m from the center of the bubble streams. In the absence of injected CO 2 , pH eddy covariance quantified the proton flux due to naturally-occurring benthic organic matter mineralization (equivalent to a dissolved inorganic carbon flux of 7.6 ± 3.3 mmol m (cid:0) 2 d (cid:0) 1 , s.e., n = 33). At the lowest injection rate, the proton flux due to CO 2 dissolution was 20-fold greater than this. To accurately quantify emission, the kinetics of the carbonate system had to be accounted for. At the peak injection rate, 73 ± 13% (s.d.) of the injected CO 2 was emitted, but when kinetics were neglected, the calculated CO 2 emission was one-fifth of this. Our results demonstrate that geochemical techniques can detect and quantify very small seafloor sources of CO 2 and attribute them to natural or abiotic origins.


Introduction
To achieve net carbon neutrality by 2050, as many governments have pledged, technology will be needed to offset CO 2 emissions during a transition to greener energy and industrial production.Carbon capture and storage (CCS) can help meet this need (IPCC, 2019).In 2019, 25 Mt of CO 2 were stored in terrestrial and offshore reservoirs (GCCSI, 2019).This is an increase over prior years but remains small compared to global CO 2 emissions of 42 Gt that year.The primary obstacles that CCS must overcome for more widespread use are insufficient public support and cost.It can be argued that public support is more important.Government subsidies for low-carbon energy such as offshore wind are similar to the subsidies that would be needed for a private CCS industry to develop (Bui et al., 2018).Highlighting the importance of public support, a terrestrial CCS project in the Netherlands was canceled due to public concern for potential environmental consequences of CO 2 leakage, among other factors (Feenstra et al., 2010).To address public concern, offshore CCS may be preferred over terrestrial CCS.In offshore CCS, CO 2 would be injected into saline aquifers or depleted oil and natural gas reservoirs hundreds to thousands of meters below the seafloor.Offshore carbon storage capacity is sufficient for mitigating emissions (e.g., Vangkilde-Pedersen et al., 2009), and offshore CCS has successfully been demonstrated at Sleipner and Snøhvit, Norway, and Tomakomai, Japan (Furre et al., 2017;Sawada et al., 2018).Nevertheless, there is public concern for the resilience of marine ecosystems to environmental damage that could be caused by offshore CCS operations (Schumann et al., 2014).Effective offshore monitoring technology can highlight issues and allow the development of strategies to mitigate risk.
Our research was conducted as part of a larger study with the aim of examining the effectiveness of chemical, acoustic, and optical sensors at detecting and quantifying a controlled release of CO 2 at the seafloor at very small rates (5.7 to 143 kg d − 1 ).For context on the magnitude of these rates, the Peterhead CCS project was designed to inject 1 million tons of CO 2 per year into the depleted Goldeneye hydrocarbon gas field at the experimental release site (Dean and Tucker, 2017).If 0.01% of this CO 2 leaked back to the atmosphere, the emission would be 270 kg d − 1 .This is the rate of CO 2 emission by a fuel-efficient car (producing 120 g CO 2 km − 1 ) travelling at ~100 km h − 1 .The goal of quantifying smaller leaks than this is to demonstrate a technical competence that can be used to ensure the integrity of storage.The rates of release are physically realistic: they are in the range of leakages that can occur along the outside of wells due to the fracturing of sediments caused while drilling (Vielstädte et al., 2019).Fortunately, the adverse environmental consequences of leaks at this low rate are minimal, with a footprint of potential harmful effects (∆pH > 0.1 unit) on the order of a hundred square meters or less (Blackford et al., 2020;Vielstädte et al., 2019).The adverse environmental consequences are also transient, with recovery occurring within a few weeks (Blackford et al., 2014).Therefore, if these leaks can be detected, then leakage into overlying water can be detected at a lower rate of emission than that which causes significant adverse environmental consequences.
A substantial challenge for the identification of a small leak is the risk of false positives, i.e., anomalously high environmental CO 2 concentrations that are attributed to CCS, but are nevertheless caused by naturally-occurring organic matter mineralization.False positives have brought accusations of leakage against a terrestrial CCS site (Beaubien et al., 2013;Romanak et al., 2013) and a temporary project suspension for offshore CCS (CSLF, 2017).Both of these events eroded public trust.Offshore CCS monitoring will face an enhanced risk of false positives due to the incentive to identify small leaks, which requires the ability to confidently differentiate a small signal from a time-varying background.Therefore, detection of CO 2 emission is only one of the tasks that will be required for monitoring.A robust CCS seafloor monitoring program would include phases of detection, attribution, quantification, and environmental assessment (Blackford et al., 2015).This program would be based on complementary technologies.To detect small leaks of CO 2 at the seafloor, active acoustics (i.e., multibeam sonar) would likely be used.Sonar detected a 29 kg d − 1 experimental CO 2 release at 110 m depth in a prior study (Dean et al., 2020) and is used for monitoring the seafloor at the Sleipner site (Linke et al., 2014).However, acoustic techniques have limitations.CO 2 bubbles, methane bubbles, and the swim bladders of fish all generate similar acoustic returns.In addition, gas concentration within bubbles cannot be resolved (e.g., Long et al., 2020), and backscatter does not correspond to bubble size, complicating quantification of emission.Therefore, acoustic techniques would likely be aided by other, complementary techniques, including optical sensors (Delwiche andHemond, 2017a, 2017b), passive acoustics (Leighton and White, 2012), and geochemical leak detection (Blackford et al., 2015).pH eddy covariance is capable of quantifying naturally-occurring benthic biological CO 2 production (Long et al., 2015).Therefore, we expect it to be highly sensitive to benthic CO 2 emission.Turbulence is the dominant mechanism responsible for the vertical exchange of mass (such as CO 2 ) with overlying air or water in terrestrial and aquatic environments.Eddy covariance quantifies the vertical flux of a solute from the covariance of high frequency (e.g., > 1 Hz), turbulent fluctuations in vertical velocity and solute concentration (Berg et al., 2003).pH and O 2 concentrations can be determined at a high-enough frequency for the technique.The dissolution of CO 2 (g) in seawater produces hydrogen ions according to Zeebe and Wolf-Gladrow (2001) as follows Because pH = -log 10 [H + ], pH eddy covariance would detect a source of CO 2 at the seafloor.
The primary goal of this study was to detect a source of CO 2 at the seafloor during a controlled, experimental release of CO 2 (g).A second goal was to quantify CO 2 released to the water column.To detect a CO 2 source (i.e., a bubble stream), we developed a pH eddy covariance system to quantify the vertical pH flux through the benthic boundary layer (the lowermost portion of the water column where current velocity is affected by bed friction).However, because CO 2 emission from the seafloor was a point source, and not spatially uniform during the experiment, eddy covariance measurements did not quantify the total CO 2 emitted.Instead, we quantified CO 2 emitted from the seafloor by calculating the advection of enriched dissolved inorganic carbon (DIC) in the plume of aqueous CO 2 generated by bubble dissolution.The mass transport of enriched DIC in the plume was calculated from current velocity, direction, and pH measured by eddy covariance instruments.Supporting measurements of high-accuracy pH and alkalinity were made by lab-on-chip sensors.The vertical distribution of DIC in the plume and its dispersion downstream were predicted by complementary numerical models.The result was adjusted to fit the observed pH at two heights above the seafloor determined by a lab-on-chip pH sensor.

Study site and the experimental release of CO 2
To provide an opportunity to demonstrate a technical competence at ensuring the integrity of offshore CO 2 storage, a controlled release of 675 kg of CO 2 (g) was conducted in the North Sea.The injection site (57 • 59.574 ′ N, 0 • 22.460 ′ W) was a candidate CCS site, a depleted hydrocarbon gas reservoir located 140 km northeast of Aberdeen, Scotland.The experiment was conducted in May of 2019 with support from two research vessels, the RRS James Cook (cruise 180) and the R/V Poseidon (cruise 534).Details of the site and the experimental design are presented in Flohr et al. (2021b).For the injection, a curved, 9-m pipe was inserted into surficial sediments to position a diffuser tip 3 m beneath the sediment surface at a distance of 7 m from the pipe inlet, which remained above the sediment surface.A remotely operated vehicle (ROV) connected the pipe inlet to CO 2 (g) delivered from a pressurized container at the seafloor.To minimize hydrodynamic interference, the CO 2 (g) container was positioned 80 m off-axis from the predominant (north-south) flow direction.CO 2 monitoring equipment including the eddy covariance instrument frame were positioned at the release site by the ROV.
The sediments overlying the pipe were sandy muds, with layers of muddy sands with a porosity of 50%.Sediments were predominantly quartz with minor calcite and clay mineral fractions (Lichtschlag et al., 2021).CO 2 was injected through the pipe continuously from 11 May to 22 May 2019.The rate of gas injection was increased step-wise, with rates of 5.7, 14.3, 28.5, 85.5, and 143 kg CO 2 d − 1 .Eddy covariance instruments were positioned 4 m to the south of the buried diffuser tip, the expected initial point of emission.Other monitoring instruments were positioned primarily to the east and west of the point of release to minimize hydrodynamic interference.

Eddy covariance measurements
Eddy covariance flux was calculated according to Berg et al. (2003) as where u z is the vertical water velocity, c is solute concentration, the prime symbol indicates a fluctuating component from which the mean has been subtracted, and the overbar indicates averaging.Water velocities were determined at 16 Hz in three dimensions with an acoustic Doppler velocimeter (ADV; Nortek Vector, Nortek AS).pH and dissolved oxygen were determined at the measurement volume of the ADV which was positioned 16 cm above the seafloor.Dissolved oxygen was measured using a non-stirring sensitive mini-optode (Holtappels et al., 2015) with a t 90 < 0.4 s (OXR-430-UHS, Pyroscience, GmbH).For LED-excitation of the sensor, and quantification of the resulting fluorescence, an oxygen meter (FSO2; Pyroscience, GmbH) was placed in a submersible housing with an optical port.Analog output from the meter D. Koopmans et al. was recorded on an analog input channel of the ADV.We chose 16 cm as the measurement height because we expected the greatest DIC enrichment just above the seafloor.Additionally, for eddy covariance measurements, selection of the measurement height is a trade-off between signal intensity and signal frequency (Long, 2021).Close to the seafloor upwards-moving eddies have dispersed less.This makes differences in their solute concentrations from the overlying water column stronger.The eddy covariance pH signal was determined by an Ion Sensitive Field Effect Transistor (ISFET; Microsens SA).Its operation relied on custom-made components following the principles of Long et al. (2015).To exclude the effect of light-sensitivity during calibrations and field work, a narrow opaque cylinder was designed to house the ISFET.Water was pumped into the cylinder and across the face of the sensor (Fig. 1).A sampling tube, 7 cm long and 2 mm in inner diameter, was used to draw water from the measurement volume of the ADV.ISFET function requires a reference voltage, for this we used a reference electrode with minimal stirring sensitivity based on a ceramic membrane design.The reference was integrated into the sensor flow-path (Fig 1).The signal of the ISFET was amplified ten-fold using an auto-zeroing galvanic isolation amplifier that has been previously used for eddy covariance O 2 measurement (Berg et al., 2003).For pH measurement, the amplifier was modified with a custom input stage for source-drain operation of the ISFET.The voltage response of the amplified ISFET signal to a change in pH was close to 590 mV per unit of pH.The precision of the amplified signal was smaller than 0.002 pH units.The amplifier was modified to compensate for the effect of temperature on the reported pH signal (25 to 40 mV per • C).To make this compensation, a temperature sensor (negative temperature coefficient thermistor) was mounted at the ISFET sensor face.Temperature measurements were used to compensate low frequency contributions (< 0.01 Hz) to the ISFET signal.To reduce noise generated by the isolation amplifier, a signal conditioning board was added to the ADV.This board also expanded the analog input range of the ADV to +/-5 V. To ensure that the sensor signal would remain in this range, the amplifier auto-zeroed the signal every half hour.Once per deployment we made a single-point calibration of the ISFET signal to background pH (8.04) determined in situ by a pH lab-on-chip sensor (accuracy +/-0.003units, precision 0.001 units; National Oceanography Centre, Southampton, UK).The lab-on-chip sensor was positioned on the same frame as the eddy covariance instruments (Fig 1), and measured pH at the same height as the ISFET (16 cm).
A reversible gear pump was developed to draw water from the measurement volume of the ADV past the ISFET.The pump gears were mounted in a 3D printed housing and operated by a motor that was programmed to reverse direction for one minute out of every thirty.The reversals cleared the pump tubing of accumulated debris.The pumped flow rate was 150 ml min − 1 .Continuous flow of water past the ceramic membrane of the reference sensor effectively eliminated its residual sensitivity to changes in water velocity (stirring sensitivity).The 90% response time of the ISFET (mean of the responses to pH increase and decrease) was 1.2 s.This is slower than the response time typically targeted with eddy covariance sensors, but it is rapid enough to record almost all of the turbulent frequencies that typically contribute to eddy covariance fluxes (e.g., 1 Hz to 0.001 Hz; Berg et al., 2013).We examined the fraction of the flux signal that was lost due to the relatively slow speed of this sensor by examining the frequencies at which contributions to the measured flux occurred during the experiment.We found that turbulent contributions to oxygen flux began at 1 Hz and turbulent contributions to pH flux began at 0.5 Hz (data not shown).In oxygen fluxes, the contribution of frequencies between 1 Hz and 0.5 Hz was less than 10% of the total signal.Therefore, pH flux underestimation due to the slower response time of this sensor was also likely to have been less than 10%.We found that during the first 48 h of eddy covariance measurements at the site, the pH signal recorded by two identical ISFET sensors was too noisy for accurate flux calculations (data not shown).The noise was unrelated to ROV operations.Bresnahan et al., (2014) identify causes of transient instability and inaccuracy that occur with the first use of a pH ISFET that likely contributed to the noise that we observed.They recommend pre-conditioning the sensor in seawater to alleviate many of these issues.
The eddy covariance pH and O 2 sensors were mounted at the measurement volume of an ADV on a lightweight fiberglass frame.An eddy covariance sensor that is placed too close to the ADV measurement volume can interfere with its velocity measurements (Berg et al., 2016).We examined the effect of sensor proximity to the measurement volume on the signal amplitude of the ADV in a tank of filtered water.Based on those measurements, we aligned the pH intake tube 2.6 cm from the edge of the measurement volume, and the tip of the O 2 sensor 1.8 cm from the edge of the measurement volume.
In benthic flux measurements, the footprint of the technique is defined as the minimum area of the seafloor that contributes 90% of the eddy covariance flux signal.It is a narrow ellipse whose dimensions are a function of the sediment roughness parameter (z 0 ) and the height of the measurement volume above the bed (Berg et al., 2007).We determined the friction velocity during northward flow over undisturbed sediments, according to (Stull, 2012) as where the x, z, and y subscripts refer to longitudinal, transverse, and vertical components of velocity direction.For this analysis, we doublerotated the velocity field of a few half-hour time intervals when flow was northward (not-influenced by the experimental site) and the signalto-noise ratio was high.Double-rotation aligned the primary axis of the velocity field with flow direction during the interval.This reduces the net transverse and vertical flow during each interval to zero.At a mean longitudinal velocity of 0.15 m s -1 , u * was 0.0116 m s − 1 .We then calculated z 0 using the law of the wall as where h is the measurement height above the seafloor, U x is the mean longitudinal velocity, and κ is the von Karman constant (equal to 0.41).
At U x equal to 0.15 m s − 1 (and h equal to 0.16 m) z 0 was 7.6 × 10 − 4 m.
Based on this measurement, the approximate area of the seafloor footprint was 30 m 2 (according to Berg et al., 2007).The predicted major and minor axes of the footprint ellipse were 35 m and 1.0 m, respectively.The point of maximum contribution to the flux signal was expected to be just 1.6 m upstream of the sensors.We report the proton flux calculated from pH eddy covariance measurements as a DIC flux.The amount of CO 2 (g) dissolution required to generate the observed proton flux was calculated using the software CO 2 SYS (Lewis and Wallace, 1998).For CO 2 SYS calculations, the seawater carbonate chemistry was constrained with measurements of pH and total alkalinity by spectrophotometric lab-on-chip sensors (Schaap et al., 2021).Additional sensors on the instrument frame included galvanic oxygen and glass electrode pH (RBR, Inc.Canada).These were also positioned to sense concentration at a height of 16 cm above the bed.Weight restrictions on the ROV limited batteries on the instrument frame; the instruments could only measure continuously for 60 h.To make near-continuous measurements of the DIC flux over the eleven days of the experiment, two identical landers were used.As one lander was retrieved from the seafloor, a second lander was positioned in its place.
The velocity, pH, and oxygen time series, recorded at 16 Hz by the ADV, were downsampled to 5 Hz, the frequency of optode data collection.Eddy covariance fluxes were calculated in half-hour intervals following the procedure described by Holtappels et al. (2013).The tilt of the ADV was corrected using the planar fit method by Wilczak et al. (2001).A running average with a window length of 300 s was subtracted from the time series to calculate the fluctuating vertical velocity (u ′ z ) and concentration (c ′ ).

DIC concentration at carbonate system equilibrium
To quantify CO 2 emission, we accounted for the kinetics of CO 2 equilibration.The kinetics of CO 2 equilibration are limited by the relatively slow reactions of the hydration and hydroxylation of CO 2 to HCO 3 − (Zeebe and Wolf-Gladrow, 2001).The equilibration time (t), defined as the time required for the carbonate system to reach 63% of equilibrium after a perturbation, is pH-dependent.As the equilibrium pH decreases, the equilibration time also decreases.To solve for the pH-dependent carbonate system equilibration time we used the complete equations presented in Schultz et al., (2006), to which we added two equations to account for its dependency on boron (Zeebe and Wolf-Gladrow, 2001).The equations were solved simultaneously using the Matlab® ODE solver.These calculations, including minor corrections to typos in prior work, are presented in Gros et al. (2021, Supplementary materials).At background pH (8.04) the equilibration time was 117.9 s.With the addition of 184 µmol L − 1 of DIC, the equilibrium pH was 7.45 and the equilibration time decreased to 85.5 s.
To calculate the concentration of excess DIC (i.e., elevation of DIC above background) at chemical equilibrium, we iteratively solved for the pH-dependent equilibration time.First, we estimated the excess hydrogen ion concentration at equilibrium (δ[H + ] eq. ).This calculation made use of the excess hydrogen ion concentration at time t according to kinetics presented by Zeebe and Wolf-Gladrow (2001) as where t is the sum of the time delay due to transport from the bubble streams to the instruments and any time required for measurement.For the first calculation, we assumed that τ = 117.9s.We then used the resulting δ[H + ] eq to calculate a new τ, according to the paragraph above.
We iterated this sequence of calculations until the δ[H + ] eq no longer changed between iterations.We then used CO 2 SYS (Lewis and Wallace, 1998) to calculate the enrichment of DIC above background concentration at equilibrium (δDIC eq ) from δ[H + ] eq .Alkalinity was 2319 µmol kg -1 (Esposito et al., 2021;Schaap et al., 2021).Our calculations were simplified by the observation that there was no alkalinity signal in the water column plume (Schaap et al., 2021).

Quantification of total CO 2 emitted
To determine CO 2 emission during the release experiment we relied on five sources of data.The first data source, eddy covariance instruments, provided time series measurements of pH, current velocity and current direction.The second data source, lab-on-chip pH sensors were used to calibrate eddy covariance pH (details in Section 2.2), and to quantify the mean vertical δDIC distribution in the plume of aqueous CO 2 generated by bubble dissolution at the maximum injection rate (details Section 3.3).Third, the predicted vertical distribution of δDIC was determined with complementary numerical models of CO 2 (g) dissolution and its dispersion (details in Section 2.5).Fourth, the distance to the center of the bubble streams was determined from ROV observations.Finally, the mean bubble diameter emerging from the seafloor was needed for model predictions.
The current direction at the release site was tidally driven in a counterclockwise motion.We divided the time series into half-hour intervals to represent separate segments of the plume.We calculated the amount of CO 2 emitted in each segment as where n refers to a segment of the plume, w is the segment width determined as a perimeter of an arc defined by the initial and final current directions during the interval, δDIC(z) is DIC enrichment above background concentration, at equilibrium, as a function of height above the seafloor, and U n (z) is the mean velocity magnitude (with x and y components) as a function of height above the seafloor.In detail, w was calculated as (θ 1 -θ 2 )⋅r, where θ 1 is the mean x-y angle of flow direction at the beginning of the interval, θ 2 is the angle at the end of the interval, and r is the radius of the implied circle.In this case, r is the distance to the primary bubble stream (2.6 m).This was the approximate distance to the center of bubble streams throughout the experiment.Bubble stream locations and their intensity were determined from photographs and video imagery of the seafloor collected by the ROV during the experiment (Flohr et al., 2021b).The DIC concentration at the height of the measurement volume was calculated as described in Section 2.3.To quantify the vertical distribution of DIC, and to investigate changes in DIC flux with distance to the plume, we used two complementary numerical models.Details are provided in Section 2.5.To quantify the vertical distribution of U n we applied the law of the wall to the mean velocity determined at the height of the ADV measurement volume.This calculation relied on the roughness parameter calculated in Eq. ( 4).DIC advection in the plume was then calculated as

Prediction of dDIC distribution using numerical models
Complementary numerical models were used to estimate the vertical distribution of DIC, and to determine DIC flux as a function of distance to the bubble streams.The first model was based on the two-phase plume model by Dewar et al. (2015).This model predicts the bubble rise height and rise velocity through a momentum balance, with the gas buoyancy counteracted by the drag encountered between the bubbles and the surrounding waters.The dissolution and distribution of DIC is predicted through a mass transfer correlation based on the bubble size, shape, and rise velocity, and the local seawater properties (i.e., pressure and temperature).The vertical DIC distribution was calculated for emerging bubbles with a mean diameter of 7 mm and a volume median diameter of 10.2 mm (determined from back-lit video imagery collected by a purpose-built lander positioned over the bubble stream on 15 May 2019 (Li et al., 2021).The second model calculated the effect of advection and dispersion on the resulting vertical distribution of DIC during its transport from the bubble-stream to the sensors.It was developed using Comsol Multiphysics® (Comsol Inc., Sweden).Specifically, we used a low Reynolds number k-ε turbulence model, coupled to the Transport of Dilute Species module.A detailed description of the underlying equations for turbulent transport, kinetic energy, its dissipation, and damping functions is provided in Holtappels et al. (2013).Briefly, the model we constructed was 12 m long in the direction of flow ( x

Leak detection with pH eddy covariance
CO 2 release commenced on 11 May 2019.As the ROV returned to the release site after opening the valves of the CO 2 container, it observed the first three CO 2 bubble streams.The most active bubble stream was located 2.6 m to the north of the eddy covariance lander.As the injection rate increased, the number of bubble streams and their areal distribution increased.Bubble streams primarily formed along a 2.5 m-wide eastwest axis centered on the original active one.An exception was a bubble stream that emerged just one meter to the north of the eddy covariance lander at the end of day 6.Its flow diminished at the peak injection rate on day 9.An image of the eddy covariance lander position relative to the dominant bubble streams on the final day of the experiment is provided (Fig. 2).By the final day, there were nine active bubble streams, two intermittent streams, and ten sites of formerly active bubble streams.Thus, there was high variability in the locations of bubble stream emission during the experiment.
The dissolution of CO 2 (g) caused a vertical, turbulent pH flux that was detected by eddy covariance instruments at the lowest release rate and throughout the experiment (DIC flux; Figs. 3, 4).The elevated DIC flux was detected on the southward portion of tidally-oscillating flow, when the instruments were exposed to the low-pH plume of the CO 2 bubble streams.At the lowest injection rate, the peak DIC flux observed in the plume of the bubble streams (> 0.020 kg m − 2 d − 1 ) was approximately 20x the background flux (< 0.001 kg m − 2 d − 1 ; Fig. 3).As the experimental injection rate increased, the peak DIC flux increased (Fig. 4).At the peak injection rate of 143 kg d − 1 the peak DIC flux was 0.5 kg m − 2 d − 1 .This is 500x the background flux, demonstrating that at close proximity, eddy covariance is remarkably sensitive to this introduced source of CO 2 .
Naively, one might expect to observe an eddy covariance DIC flux of the rate of injection divided by the area of emission.At the peak injection rate of 143 kg d − 1 , the bubble streams could be contained within a circle with a radius of 1.75 m.This is an area of 9.6 m 2 , giving a firstorder predicted emission of 14.9 kg m − 2 d − 1 (143 kg d − 1 divided by 9.6 m 2 ).Instead, we observed a mean DIC flux of 0.06 kg m − 2 d − 1 when the instruments were downstream of the bubble streams (Fig. 4).A combination of factors contributes to this small flux.First, we can expect that a significant fraction of the injected CO 2 was retained by sediments (see Section 4.4).Second, of the CO 2 (g) that is emitted to the water column, much of it dissolves above the measurement volume of the eddy covariance instruments and therefore remains undetected by eddy covariance (Dewar et al., 2015;Gros et al., 2021).Third, the footprint of the eddy covariance technique is long, narrow, and larger than the area over which CO 2 emission occurred.Finally, the low pH plume generated by CO 2 bubble dissolution was likely not at chemical equilibrium due to the close proximity of instruments to the bubble streams (see Section 3.4).

Biotic CO 2 (aq) production and O 2 uptake
The pH eddy covariance technique was sensitive enough to quantify the pH flux that results from biological CO 2 production in seafloor sediments at the site.For these measurements we examined pH flux during northward flow, when the eddy covariance footprint was entirely upstream of the bubble streams.We collected the highest quality fluxes (i.e., fluxes with the most consistent flux signal) on day ten of the release experiment.The mean biological DIC efflux matched the mean O 2 uptake determined at the same time (Fig. 5).Differences between biological O 2 uptake and DIC production were not statistically significant (ANOVA, α = 0.05), but both fluxes were significantly different from zero (one-sample t-tests, α = 0.05).ROV disturbance of sediments on the south side of the lander during lander deployment and retrieval may have contributed to noise in the flux signal, but not all of the variation in biotic flux was noise.The magnitude of oxygen uptake and DIC production increased with increases in velocity (e.g., Fig. 5 at day 10.3).Benthic heterotrophic respiration is ultimately limited by organic matter delivery and its quality (Arndt et al., 2013).Over short time-scales, however, respiration may be enhanced by increases in oxygen supply from overlying water.Increases in velocity can also lead to increases in solute exchange.For these reasons, tidally driven oscillations in oxygen uptake occur in cohesive sediments (e.g., Donis et al., 2016;Glud et al., 2016;Koopmans et al., 2021).Further research is needed, but our results suggest that biotic DIC fluxes are similarly enhanced over short time scales by increases in water velocity.
The high sensitivity of the eddy covariance technique results from the covariance of high frequency velocity and pH signals.The covariance is illustrated in Fig. 6.At the peak injection rate, pH decreases greater than 0.03 units occurred roughly every four seconds in the plume of aqueous CO 2 generated by the experimental injection.Decreases in pH commonly occurred with an upward vertical velocity.The converse was also common.Both cases result in a net upward flux of hydrogen ions.Thus, both cases contribute to calculated DIC flux (Eqn. (1)).The covariance of vertical velocity and pH shows that the high frequency variations are a flux signal, not noise.

DIC flux as a function of injection rate
The peak eddy covariance DIC flux in the plume of the bubble streams did not increase linearly with the injection rate (Fig. 7).At injection rates above 20 kg d − 1 , the additional CO 2 injected resulted in smaller increases in DIC flux.An exception is one observation at the highest injection rate.At this time, the peak DIC flux varied up to fourfold during the 13 h elapsed from one observation of the plume to the next.The diminished response above 20 kg d − 1 was caused by the increasing area and number of bubble streams through which CO 2 emission occurred at the higher injection rates.Only a subset of the bubble streams would contribute to measured eddy covariance fluxes at any one time.The high variability in flux at the highest injection rate was likely due to intermittent CO 2 emission by the closest bubble streams.
To quantify the effect of bubble stream proximity on DIC flux, and to examine the vertical and horizontal distribution of bubble-derived DIC, we used complementary numerical models.The two-phase model predicted that the maximum rate of CO 2 (g) dissolution occurred just above the seafloor (Fig. 8A).This is where the volume of CO 2 (g) is at a maximum.The predicted rate of CO 2 (g) dissolution diminished nearlinearly with height and, for a mean initial bubble diameter of 7 mm, was negligible at heights of greater than 4 m.The vertical profile of water velocity is inverted relative to this, with minimal velocity and minimal turbulent mixing near the bed.As a result, the advectiondispersion model predicted that DIC enrichment would be close to its maximum near the height of our instruments (0.16 m; the lower of the observation points in Fig. 8B).Lab-on-chip measurements of the vertical pH gradient at two points (0.16 m and 0.88 m) revealed a steeper mean gradient at the maximum release rate than predicted by our models.The cause may be an underestimation of the abundance of small bubbles (e. g., < 2 mm diameter), which would dissolve at a low height, enriching    δDIC near the seafloor.Underestimation of small bubble sizes, along with coalescence and breakup of bubbles at low depths, make in situ bubble measurement challenging (Sellami et al., 2015).To improve the accuracy of our quantification of the vertical DIC profile, we fit our prediction to the observed concentration gradient (Fig. 8B).The resulting gradient is similar to one predicted by the emergence of bubbles with a mean diameter of 4 mm (Schaap et al., 2021).
The advection-dispersion model predicted that our measured eddy covariance DIC flux would depend on the distance to a bubble stream (Fig. 8C).At a distance of 1 m, the flux would be twice as great as at 2.6 m.Therefore, the temporary appearance of a bubble stream close to the instruments can explain the high variance in eddy covariance DIC flux at the peak injection rate Fig. 7.This model suggests that eddy covariance instruments would need to be located close to the bubble streams (i.e., within 10 m) in order to detect the smallest release rate in this study.

Equilibrium of the carbonate system
We found a substantial discrepancy between eddy covariance pH and lab-on-chip pH in the plume of the bubble streams.At the peak injection rate, the eddy covariance pH sensor and an accompanying glass pH electrode both reported a mean pH of 8.01 during southbound currents.The lab-on-chip pH sensor, measuring at the same height and on the same frame, reported a mean pH of 7.83 (e.g., Fig. 9A).The difference was due to incomplete chemical equilibration of CO 2 .According to our calculations the equilibration time of the carbonate system was 117.9 s (at pH 8.04, 7.8 • C).At a mean water velocity in the plume of 11 cm s − 1 , and a distance of 2.6 m between the bubble streams and the instruments, the eddy covariance sensor and glass electrode would have detected the pH signal of CO 2 dissolution 24 s after the CO 2 bubble emerged from the sediment.This short time would allow for only 18% of CO 2 hydration/ hydroxylation and proton production (Eq.( 6)).In contrast to the eddy covariance sensor, the lab-on-chip sensor includes an average time-lag of 140 s between sample collection and analysis.Because of this, labon-chip pH would have been determined at 70% of chemical equilibrium.
For this study, we used pH observations made during eddy covariance pump reversals to determine δDIC at chemical equilibrium.Lab-onchip pH measurements were unavailable for this because they are the subject of a companion manuscript (Schaap et al., 2021).During standard eddy covariance measurements, inflowing water passes the ISFET sensor and travels through 6 cm 3 of tubing on route to a 550 cm 3 filter housing that traps sand before reaching the gear pump.During the one-minute-long reversal, 150 cm 3 of this water are discharged backwards past the ISFET sensor, providing a second pH measurement.This temporary storage of water brings the carbonate system closer to chemical equilibrium.The pH of discharging water in the final seconds of pump reversal, hereafter referred to as the pump reversal pH, was similar to observations of pH detected by the lab-on-chip sensor (Fig. 9A).At the peak injection rate, the mean pump reversal pH in the plume of the bubble stream was 7.86, close to the mean LOC pH of 7.83.Outside of the plume, the pump reversal pH matched the background (8.04), indicating that the reduced pH detected in the plume was due to CO 2 (g) dissolution caused by the bubble streams.We used the successive measurements (inflow and outflow) to constrain the transport time and residence time of water in the system.We found the best fit of mean residence time of water in the system (data not shown) was 105 s (the midpoint of inflow pH shown in Fig. 9B).We then used pump reversal pH to calculate δDIC at equilibrium (Fig. 9C).

Quantification of CO 2 (g) emission
Time series measurements of velocity and pH were used to calculate the transport of δDIC in narrow arc-lengths defined by half-hour measurement intervals according to Eq. ( 6).An illustration of the results of these calculations is presented in Fig. 10.The equilibrium δDIC content of the plume was highly heterogeneous in cross-section, with differences of 100 mmol δDIC L − 1 separated by tens of centimeters.Based on our calculations of chemical equilibrium, 73 ± 13% (mean ± s.d.) of CO 2 was emitted to overlying water at the peak injection rate (Fig. 11).Uncertainty in these measurements is the normalized s.d. of emission at  injection rates of 88 and 143 kg d − 1 .If we disregard chemical equilibrium and base our calculations of emission on pH measured during inflow, the calculated emission was one-fifth of this.Therefore, to accurately determine CO 2 emission close to a bubble stream, the equilibration of the carbonate system must be accounted for.

Leak detection with pH eddy covariance
Our results demonstrate that pH eddy covariance is an effective technique for detecting extremely small CO 2 leaks from the seafloor.The limit of detection in other settings would depend on distance to the leak, instrument height, and water depth, among other factors.However, we unambiguously detected CO 2 (g) bubble dissolution at the lowest injection rate of 5.7 kg d − 1 (Fig. 3).We were also able to quantify biotic CO 2 (g) emission (Fig 5).Active acoustic techniques are a highly effective technology at bubble detection, and are currently in use for monitoring of offshore CCS sites.We suggest that pH eddy covariance could be useful in a complementary monitoring role because of its high sensitivity and chemical specificity.A hypothetical application would be to discriminate methane from CO 2 bubble streams.Eddy covariance could also be used to detect dissolved CO 2 emission where bubbles are absent, such as at great Among the restrictions on the use of eddy covariance in these roles is a requirement that the instruments may need to be within 10 m of the source to detect small leaks (e.g., 5.7 kg d − 1 ; Fig. 8).The point of maximum contribution to eddy covariance flux during typical deployments is within 1-2 m of the instruments (Berg et al., 2007;Rodil et al., 2019).In this respect, while eddy covariance appears to be more sensitive at leak detection than other geochemical techniques, it may not substantially extend the expected distance to which a CO 2 bubble stream could be detected by geochemical techniques.This limitation of proximity could be overcome by mounting eddy covariance instruments to an autonomous underwater vehicle.Aquatic eddy covariance has previously been determined from a moving platform (Berg et al., 2020;Flügge et al., 2016;Long and Nicholson, 2018).Therefore, the method is feasible.Alternately, pH eddy covariance could be used to monitor a location at risk of leakage, such as at a well-head.

Importance of equilibration of carbonate chemistry
The difference in DIC emission calculated at the pH of inflow and at chemical equilibrium (Fig. 11) demonstrates the importance of accounting for chemical equilibrium when measuring CO 2 emission close to a bubble stream.At the temperature of the experimental release site, and at a mean water velocity of 11 cm s − 1 , the carbonate system would be at only 79% of equilibrium at a distance of 20 m (calculated following Eq.( 5), where τ is 117.9 s).Therefore, to a measurement distance of tens of meters, the effect can be significant.The time to equilibration diminishes rapidly as temperature increases.At 16 • C, for example, τ is less than 50 s.Therefore, at higher temperatures this effect would be less of a concern.

Biotic CO 2 (aq) production and O 2 uptake
During measurements of biotic CO 2 production using the pH eddy covariance technique, biogeochemical processes other than CO 2 (g) dissolution may contribute to proton flux.In broad terms, sediment pH is driven by aerobic and anaerobic organic matter mineralization as well as calcium carbonate dissolution and precipitation (reviewed by Burdige, 2006).A simplifying observation is that whether the perturbance of the carbonate system is due to production of protons or production of CO 2 , the equilibration of the carbonate system drives an equivalence between them, such that a proton flux represents a net CO 2 flux (see eqn. ( 1)).
We suggest that disequilibrium of carbonate chemistry is not an issue for calculations of benthic, biological (biotic) CO 2 production.A fraction of benthic biotic CO 2 production may escape detection in pH eddy covariance measurements because it has not yet hydrated/hyroxylated in water, i.e., it was produced within 2-3 min of it passing across the face of the pH sensor.A complete analysis is beyond the scope of this work,  but we suggest that the fraction that escapes detection is too small to alter our calculated biotic CO 2 emission.Diffusive transport times for solutes through sediments and the diffusive boundary layer can be calculated according to (Jørgensen and Revsbech, 1985) as where z is the diffusion distance, and D is the molecular diffusivity (10 − 5 cm 2 s − 1 for CO 2 at 8 • C).For CO 2 produced just below the sediment surface the diffusion distance might be just 400 µm (e.g., 100 µm of pore space and 300 µm of diffusive boundary layer at the sediment-water interface).In this case, the diffusive transport time would be 141 s.If this occurred in sediments close to the eddy covariance instruments the turbulent transport time could be as short as 10 s.The diffusive and turbulent transport times would then add to 151 s.This minimal transport time allows 72% of CO 2 to hydrate/hydroxylate prior to measurement (Eq.( 5)).The vast majority of benthic, biotic CO 2 production occurs at greater depths than this (e.g., Jørgensen, 1982), and diffusive transport time increases with z 2 (Eq.( 9)).Therefore, in most cases, benthic, biotically produced CO 2 would equilibrate in water prior to passing across an eddy covariance sensor.Biotic CO 2 production in sediments, calculated from H + flux, was similar to O 2 uptake (Fig. 5).Generally, one would expect CO 2 production to match O 2 uptake during organic matter mineralization in marine sediments.This is because processes that consume oxygen without producing CO 2 tend to be offset by processes that produce CO 2 without consuming oxygen.Specifically, the oxidation of non-carbon elements in organic matter consumes additional oxygen.This oxygen consumption is often offset by denitrification and iron pyrite formation in which N 2 and FeS 2 are the ultimate electron acceptors (Glud, 2008).CaCO 3 production or dissolution in sediments can substantially alter this ratio.CaCO 3 dissolution is a source of CO 2 but a sink for H + .For our H + flux measurements, net CaCO 3 dissolution in sediments would diminish our calculated CO 2 emission.In the silicate sediments at the experimental release site, however, net CaCO 3 production or dissolution was minimal (Dale et al., 2021).

Quantification of CO 2 (g) emission
In terrestrial CCS, eddy covariance has been used to quantify experimental CO 2 release (Lewicki and Hilley, 2009;Lewicki et al., 2010).Those studies relied on a more advanced approach for quantification, a least squares inversion of footprint functions.In our study, we were fortunate that the regular, counterclockwise progression of current directions, and the rapid dissolution of CO 2 into water allowed us to constrain plume dimensions from measurements taken at a single point.Among the limitations of our approach is that the distance to the bubble stream determines at a ratio of 1:1 the transverse dimension (w n ) of the plume, and therefore its magnitude.Our measurement of CO 2 emission at the peak CO 2 (g) injection rate (73 ± 13%) indicates that 27% of the injected CO 2 (g) was retained in sediments at that time.A total of 675 kg of CO 2 were injected during the experiment.If we apply the 27% loss rate to the whole experiment, then 183 kg were retained in sediments during this experiment.Given that the depth of the diffuser was only 3 m, the results suggest that sediments were a very effective sink for CO 2 (g).Companion studies suggest that much of the lost CO 2 was dissolved in porewater.Evidence for this includes a 0.5 • C temperature increase at 10 cm depth in the vent area due to CO 2 -induced carbonate and silicate mineral dissolution (de Beer et al., 2021).Alkalinity accumulated in sediments as a result of this dissolution (Lichtschlag et al., 2021) showed evidence of gas pockets within the sediments (Roche et al., 2021).These results are consistent with substantial CO 2 loss to sediments in a similar experiment during which CO 2 was injected at greater depths at a nearshore marine site (Lichtschlag et al., 2015).
Our results compare favorably with measurements of CO 2 release to the water column made by other techniques during this experiment.In closely related work Schaap et al. (2021) used lab-on-chip sensors to quantify the vertical pH gradient every 20 min in the plume of aqueous CO 2 generated by the bubble streams.They found that after accounting for equilibration, 61 ± 10% of injected CO 2 was released to the water column at the peak injection rate.In another geochemical approach, Gros et al. (2021) used pCO 2 sensors mounted to a CTD to survey the plume of aqueous CO 2 .They found that 64% of injected CO 2 was released to the water column at the peak injection rate.(Flohr et al., 2021a) quantified release to the water column by gas collection over bubble streams.They found that 38 ± 12% and 48 ± 7% of injected CO 2 were released to the water column during successive days at the peak injection rate.Finally, hydrophones were used to quantify CO 2 bubble emission from the "pop" sound they make as surface tension encloses a bubble on its emergence from the sediment (Li et al., 2021).They found that 22 ± 62% of injected CO 2 was released to the water column at the peak injection rate.

Anomaly attribution
Beyond leak detection and quantification, eddy covariance may be useful for attribution of leak-like pH signals in a seafloor CCS monitoring program.Geochemical leak detection has focused on identification of variations in pH with flow direction that can be indicative of a leak (e.g., Blackford et al., 2017;Oleynik et al., 2018).However, natural processes can also cause variations in pH with flow direction.Prior to the start of CO 2 injection in this study, tidally-oscillating pH variations of 0.01 unit were observed by lab-on-chip pH sensors at the release site (Schaap et al., 2021).These oscillations were caused by velocity-induced redistribution of the equilibrium pH profile in the benthic boundary layer.The effect of transient velocity on equilibrium dissolved oxygen profiles has been previously described (Holtappels et al., 2013).Stable stratification could also expose stationary sensors to variations in pH with flow direction.Similar oscillations in pH may also be caused by internal waves.In any of these cases, eddy covariance could improve attribution of anomalies like these by quantifying solute flux instead of solute concentration.Additionally, the abiotic or biotic origin of a CO 2 source can be determined by comparing eddy covariance oxygen uptake to DIC efflux.Where DIC efflux exceeds uptake, the source is abiotic.This work would complement other approaches that resolve the concentration changes of biologically-linked elements to identify anomalies caused by leaks (e.g., oxygen consumption; Botnen et al., 2015;Omar et al., 2021).Improved models of the physics, chemistry, and biological element cycling in seawater overlying offshore CCS sites would improve our understanding of these dynamics.Here, eddy covariance can also contribute by improving quantification of biotic fluxes and fluxes across pycnoclines (e.g., Kreling et al., 2014;Rovelli et al., 2016), improving the accuracy of models.

Conclusion
Eddy covariance may be useful for offshore CCS, including environmental monitoring, leak detection, attribution of leak-like signals, and quantification of CO 2 emission.The most remarkable attribute may be that because natural CO 2 production can be quantified with the technique, we can be confident that eddy covariance can detect small leaks (e.g., 5.7 kg d − 1 in this study).The covariance of vertical velocity and pH makes the eddy covariance technique, in effect, an amplifier for the detection of CO 2 (g) dissolution at the seafloor.The combination of velocity and pH sensors can also be used to quantify the total CO 2 (g) emission from the seafloor using the plume advection approach that we demonstrated.For quantification, the distance to the source and an estimate of the distribution of CO 2 (g) dissolution with height would be needed.To be useful at leak-detection in CCS monitoring, the instruments could be placed to monitor at-risk locations such as a well head, or positioned on mobile platforms.D. Koopmans et al.

Fig. 1 .
Fig. 1.Design and deployment of the pH eddy covariance system.A) Illustration of the pH ISFET housing, integrated AgCl reference, and its ceramic membrane (yellow) in the path of pumped flow.B) Arrangement of the intake of the pH sensor and the tip of the O 2 sensor outside the measurement volume of the acoustic Doppler velocimeter.C) Eddy covariance hardware and lab-onchip sensors mounted to a fiberglass frame for deployment.

Fig. 2 .
Fig. 2.The release experiment site.Top) Eddy covariance instruments and lab-on-chip sensors positioned a few meters to the south of CO 2 bubble streams on day 10 of the release experiment (21 May 2019, injection rate of 143 kg d − 1 ).Bottom) Schematic of instrument arrangement at the site modified fromFlohr et al., (2021).The EC/gradients lander (red) is the eddy covariance and lab-on-chip lander.

Fig. 3 .
Fig. 3. Detailed view of the detection of CO 2 (g) dissolution (DIC flux) at the lowest injection rate (5.7 kg d − 1 ).Top) Current magnitude and direction with respect to north.Center) High frequency and time-averaged eddy covariance pH measurements.Bottom) DIC flux calculated from eddy covariance pH flux.Throughout this manuscript DIC flux is reported in kg of CO 2 for comparison with CO 2 injected.

Fig. 4 .
Fig. 4. DIC flux calculated from eddy covariance pH flux during the release experiment.Top) Current magnitude and direction with respect to north.Bottom) DIC flux and the experimental CO 2 (g) injection rate.The gap in current and DIC flux at day 9 was due to biofouling.Note that fluxes on days 2 and 3 are also presented in Fig. 3.

Fig. 6 .
Fig. 6.Example of the covariance of turbulent vertical velocity and [H + ] that gives eddy covariance its high sensitivity.Top) The turbulent, fluctuating component of vertical velocity.Center) The turbulent, fluctuating component of [H + ].Bottom) Cumulative flux contributions occur almost continuously in the plume of the bubble stream (day 10).

Fig. 7 .
Fig. 7.The peak eddy covariance DIC flux (mean ± s.d.) in the plumes of the bubble streams.The three obervations of peak flux at the peak injection rate are shown (open circles).Biotic CO 2 production is represented at an injection rate of zero.

Fig. 8 .
Fig. 8. Model prediction of A) the dissolution of CO 2 (g) by bubbles during this experiment, and B) the resulting vertical gradient in excess DIC at the instrument location fit to lab-on-chip pH observations.C) Prediction of DIC flux with distance downstream of a bubble stream.The modeled measuring height was 0.16 m.Flux is less predictable at the bubble stream (dotted line).

Fig. 9 .
Fig. 9. Disequilibrium of carbonate chemistry in the plume of the bubble streams at the peak CO 2 injection rate.A) pH determined by sensors that briefly retained a water sample prior to analysis observed a lower pH.B) Pumpreversal pH measurements on the same parcel of water during inflow (90 s running average) and at the end of the reversal at points L and R in panel A. C) δDIC at inflow pH, pump-reversal pH, and at chemical equilibrium at points L and R in panel A.

Fig. 10 .
Fig. 10.Example calculation of plume δDIC transport.A) Time series of current magnitude, direction, and pH determined by the eddy sensor during inflow (light gray, with 5-minute running average in dark gray), during pump-reversal, and at chemical equilibrium.B) The resulting δDIC distribution and transport in cross-section.