The ocean acidification problem is becoming more and more serious with the continuous increase of atmospheric CO₂ from fossil fuel burning. Ocean acidification can be defined as the increase of dissolved CO₂ (CO₂(aq)) and consequent decrease of pH in seawater, with increases of dissolved inorganic carbon (DIC) but little variations in total alkalinity (Gattuso and Hansson, 2011). In the last two decades, thousands of studies have been carried out to study the ocean acidification effects on different marine organisms which have been reviewed and synthesized (e.g., Hoegh-Guldberg et al., 2007; Riebesell and Tortell, 2011; Meyer and Riebesell, 2015; Lemasson et al., 2017). These studies have shown that ocean acidification has complex effects on marine calcifiers (Figuerola et al., 2021), non-calcifying marine life (Hurd et al., 2019), and therefore profound impact on marine ecosystems and ocean carbon cycles (Mostofa et al., 2016). Over the past decade, more studies have employed isotopic methods in laboratory cultures to trace the ratio of stable isotopes. The variations of stable isotopes could reveal important physiological responses to ocean acidification beyond, for instance, growth rate, cell size, or elemental stoichiometry. Moreover, predicting these isotopic variations is crucial to calibrate new proxies for reconstructing the atmospheric CO₂ concentration in geological history (Hopkinson et al., 2011; Wilkes et al., 2017; Nishida et al., 2020; Phelps et al., 2021; Remize et al., 2021).

Laboratory culture is a key method to study the physiological effect of ocean acidification on different marine life. There are multiple methods to achieve the target culture media CO₂(aq) and carbonate chemistry depending on the objectives of the study. The principal methods are (1) manipulating pH by adding acid/base, (2) manipulating DIC through addition of HCO₃ ^- or CO₃ ^2- and (3) bubbling (or aeration) a gas of desired pCO₂ concentration (Gattuso et al., 2010). The method of bubbling cultures with CO₂ requires gas mixing to obtain the desired CO₂ level and humidification to prevent evaporation from culture media. Adding acid/base removes the mechanical stress upon cells from bubbles and benefits from relatively easy operations; however, it could cause shifts of alkalinity and does not reproduce the actual mechanism by which DIC concentration increases in the ongoing Ocean Acidification. Both bubbling and acid-base manipulation can well simulate the CO₂ increase and pH decrease effects in laboratory culture, but the CO₂ bubbling method has been preferred by some studies because it alone can perfectly replicate the current ocean acidification caused by anthropogenic CO₂ without changing the seawater total alkalinity.

For the CO₂ bubbling method, the guidebook by Riebesell et al. (2011), covering the methods of laboratory culture for ocean acidification research, highlighted the importance of pre-equilibrating the culture media to the required CO₂ concentration by aerating it ‘for a few days’. Considering the wide range of culture vessel shapes and volumes among experiments, our question is how long culture media should be bubbled in order to reach an ordinary chemical and isotopic balance. Some published works mentioned the pre-bubbling durations, for example, the seawater was pre-bubbled for 2 days (Iglesias-Rodriguez et al., 2008), while most publications did not fully describe their methods. Moreover, the isotopic equilibration times are usually much longer than the ordinary chemical equilibration times, because, to reach isotopic equilibrium, each ion and molecule should be fully exchanged and come to equilibrium with other ions and molecules (Mills and Urey, 1940). For the works focusing on organic or carbonate carbon isotope fractionations under different CO₂ levels, culture media with out of equilibrium or dynamic carbon isotope ratio of DIC could complicate or even preclude the interpretation of stable isotope fractionation signatures.

In this study, we provide a thorough characterization of the isotopic equilibration process in CO₂ bubbling experiments and the factors that influence the carbon isotopic equilibration time, in order to clearly document the approaches needed to accurately infer carbon isotopic fractionations in experiments with bubbling. First, we compose numerical models to simulate chemical and isotopic equilibration during bubbling processes in two different systems and present the effects of seawater volume, gas exchange rate, and difference between the isotopic composition of the bubbled CO₂ and that of the media’s DIC on the equilibration time. Secondly, we complete a series of bubbling experiments in a photobioreactor to test the performance of the model simulation. Finally, we evaluate the expected consequences of equilibration time in typical experimental bubbling setups for which carbon isotopic ratios of cultured biomass or biominerals have been reported. With the guide of this study, future works can trace the isotopic fingerprint of ocean acidification on marine biomass better.

We have conducted experiments in commercial photobioreactors of 1 L and 3 L (PBR FMT 150, Photon Systems International) designed for continuous culture. The aeration system allowed gas to first enter the bottom of a bottle with fresh media (‘bottle’ in Figure 1A ), where gas humidification and the first exchange of gases occurred. The gas subsequently flowed out towards the photobioreactor where a sparging tube dispersed the bubbles, exchanging gases a second time. Finally, gas flowed to the waste bottle, and from there out of the system ( Figure 1A ). The photobioreactor compartments were monitored without inoculated cells in batch mode in order to assess the dynamics of the bubbling process itself. The two different photobioreactor sizes and CO₂ concentrations employed in this study can be found in Table 1 .

For the different CO₂ concentration treatments, two compressed gases, pure CO₂ with a δ¹³C of -2.8‰ (Vienna Pee Dee Belemnite, VPDB) and CO₂-free synthetic air (Air Liquide), were mixed with a Gas Mixing System (GMS-150, Photon Systems Instruments). GMS output flow and pCO₂ composition were further verified with a flowmeter and a cavity ringdown spectrometer isotopic and gas concentration analyzer (G2131-i, Picarro, Inc. USA). Measurements of the Picarro CO₂ analyzer were corrected with CO₂ mixtures with certified concentrations and isotopic composition (Air Liquide).

Photobioreactors and bottles were filled with K/2 media without Tris buffer which was made by artificial seawater (ASW) with nutrient and trace metals (Keller et al., 1987). For the convenience of description, the terms of media, ASW and seawater will be treated equally to the K/2 media in the following without specific mention. The detailed recipe is described in Appendix A . Prior to filtration, 10 L artificial seawater (ASW) was supplemented with 20.48 mL 2M Na₂CO₃ and 2.90 ml 37% HCl to raise alkalinity above 4 mmol kg^-1 seawater. This operation of increasing alkalinity is specialized to maintain the growth pH in high CO₂ treatment. The salinity of ASW would increase slightly from 34.97 to 35.53 g kg^-1 and this increment in salinity could be ignored during the carbon system calculations.

Concentration of DIC in seawater and concentration and carbon isotope ratio of CO₂ in headspace were monitored with an Apollo SciTech DIC-C13 Analyzer coupled to the Picarro CO₂ analyzer using in-house NaHCO₃ standards dissolved in deionized water at different known concentrations and δ¹³C values from -4.66 to -7.94‰. δ¹³C-DIC in media were measured with a Gas Bench II with an autosampler (CTC Analytics AG, Switzerland) coupled to ConFlow IV Interface and a Delta V Plus mass spectrometer (Thermo Fischer Scientific). The system and abovementioned in-house standards were calibrated using international standards NBS 18 (-5.014‰) and NBS 19 (+1.95‰). The analytical error for CO₂(g) concentration is <20 ppm and that for DIC concentration and δ¹³C is <10 μM and 0.1‰, respectively.

Our initial δ¹³C of un-bubbled DIC (at t₀) is -6.1 ± 0.2‰. Before and after the start of bubbling at a flow of 200 ± 20 mL min^-1, both headspace and seawater media of the upstream bottle and the photobioreactor were sampled by a 50-mL syringe through one-way sampling ports. The sampling time in each experiment can be found in Tables 2 , 3 . To measure headspace CO₂(g) that had been humidified and exchanged with bottle media, gas flow was directed into a syringe and 50 mL of gas were injected into the Picarro CO₂ analyzer.

To measure seawater DIC, pH and δ¹³C_DIC, 35 mL seawater were syringed out as depicted in Figure 1A . The first 5-10 mL out of 35 mL were routinely discarded to avoid mixing effects with dead volumes in the tubing. One mL was injected into He-flushed glass vials containing H₃PO₄ for the carbon isotope ratios measurements in Gas Bench. About 12 mL were injected into the glass vials without headspace, for DIC concentration measurements using Apollo analyzer, and 3.5 mL was consumed in each duplicate. There was no gas exchange between air and samples during the sampling of DIC and CO₂. The remaining seawater in syringe (~15-20 mL) was then injected into a 50 mL centrifuge tube for pH measurement. The tube was covered with Parafilm M sealing film to reduce the CO₂ exchanging between seawater and air. The pH of seawater was measured by a Mettler Toledo LE410 pH-probe calibrated with NBS standards (Mettler Toledo) resulting in an accuracy of ±0.01. Here it should be noted that high ionic strength calibration standards, such as standards in total scale, would be optimal for pH measurement in seawater (Kadis and Leito, 2010).

There are three important processes in DIC carbon isotope evolution simulations, (1) CO₂ exchanging between gas (CO₂(g)) and seawater (CO₂(aq)), (2) DIC inter-reactions and (3) isotopic fractionation during the DIC reactions, which will be introduced separately in the following sections. Beside these three main processes, the sampling of DIC and gas in headspace can also play a minor role in DIC isotope evolution by decreasing the total amount of DIC and accelerating isotopic equilibrium. Thus, the decreasing of DIC volume and the loses of CO₂(g) in headspace during sampling are also considered in our model. As described in the last section, the CO₂ coming from the Gas Mixing System first goes into the bubbling in bottle, exchanging with DIC in bottle. Then CO₂(g) goes out of the seawater in bottle into the bottle headspace. After that, CO₂(g) goes into bubbles in photobioreactor exchanging with DIC in photobioreactor. However, in our model, bubbles in bottle and photobioreactor are combined with bottle headspace to reduce the calculation amount. Thereby, in practice, the simulated CO₂(g) goes into headspace directly after flowing out of Gas Mixing System, and exchanges with DIC in bottle and photobioreactor together ( Figure 1B ). With these simplifications, there are only two degrees of freedom in our model: CO₂ exchange rate constants (k_E) in bottle and photobioreactor. Using a given combination of k_E, the forward model runs ordinary differential equations (ODEs) toward steady state using the Matlab function ‘ode15s’, with seawater and CO₂(g) composition in bottle, photobioreactor and bottle headspace as initial conditions. The notations and equations of the model are described in detail in the Appendix B, C , respectively. Fitting processes were carried out to estimate the exchange rate constants and gas flux. These processes were achieved by minimizing the difference between simulated carbon isotope ratios and measured values via the Matlab function ‘fmincon’.

In this study, we carried out two experiments to estimate the CO₂ exchange rate constants between gas and seawater. The fitting results of CO₂ exchange rate constant (k_E) are 8.71×10^-5, 4.57×10^-5, 5.32×10^-5 and 3.36×10^-5 mol s^-1 atm^-1 for large system bottle, large system photobioreactor, small system bottle and small system photobioreactor, respectively.

The δ¹³C_DIC before bubbling are around -6.1‰ (-6.54~-5.65‰). With the onset of bubbling, δ¹³C_DIC responded logarithmically, increasing fastest during the first hours and slowing the rate of increase in the following days. The δ¹³C_DIC in both experiments did not increase further after reaching values around 5.5 ‰, about 8.3 ‰ higher than the CO₂(g), which well fitted our prediction in the last section. The δ¹³C_DIC reached equilibrium with CO₂(g) at 6 days in low CO₂ experiment with pCO₂ = 470 ppm, while in the other experiment, the isotopic equilibrium was achieved at 5 days after bubbling. In our simulations, the carbon exchange rate between CO₂(aq) and HCO₃ ^- is more than two orders of magnitude higher than the rate between CO₂(g) and CO₂(aq). Therefore, carbon isotope ratios of CO₂(aq) (δ¹³C_CO2(aq)) are almost parallel with δ¹³C_DIC (dashed lines in Figure 3 ).

Compared to the continually increasing δ¹³C_DIC and δ¹³C_CO2(aq), the carbon isotope ratio of CO₂ gas (δ¹³C_CO2(g)) in bottle headspace interestingly showed more variations (blue dots in Figure 3 ). The initial value of δ¹³C_CO2(g) was around -15‰, which is the atmosphere CO₂ carbon isotope ratio in the poorly ventilated laboratory. There were sharp increases in δ¹³C_CO2(g) from -15‰ to around -6‰, immediately after bubbling (-4.31‰ in high pCO₂ experiment and -7.32‰ in low pCO₂ experiment, only five minutes after bubbling). This was caused by the CO₂ in bottle headspace being rapidly replaced by the new CO₂ coming from the Gas Mixing System, which has a carbon isotope fingerprint of -2.8‰. With a fixed gas flux, this kind of rapid increase in δ¹³C_CO2(g) was more significant in high CO₂ concentration experiment ( Figure 3B ). The rapid increase of δ¹³C_CO2(g) was then followed by a decrease of δ¹³C_CO2(g), which was caused by CO₂ exchanging between gas in headspace and DIC in seawater. In the large system, the CO₂ exchange rate is about 640% higher than the rate in small system, due to the higher pCO₂ and higher k_E. Therefore, the ¹³C in CO₂(g) went into DIC in seawater faster in the larger system, resulting in an about 1.2‰ decline in δ¹³C_CO2(g) and also faster increases in δ¹³C_DIC ( Figure 3B ). This complex pattern of δ¹³C_CO2(g) was well simulated in our model (blue lines in Figure 3 ), though the simulation results are a bit lower value than measurements in high CO₂ experiment. This could be caused by the combination of the bubbles in bottle and photobioreactor with the headspace in our model, resulting in a more significant decline in δ¹³C_CO2(g) when CO₂(g) begins to exchange with DIC.

The original contributions presented in the study are included in the article/ Supplementary Material . Further inquiries can be directed to the corresponding author.

IT-R and HZ carried out the bubbling experiments. IT-R and PA measured carbon isotope with help from MJ in data calibration. HZ developed the numerical model. HZ and HS wrote the paper with input from other authors. All authors contributed to the article and approved the submitted version.