Contrasting responses of the ocean's oxygen minimum zones to artificial re-oxygenation

We used the Nucleus for European Modelling of the Ocean (NEMO) platform (version 4.0). The NEMO platform comprises an ocean dynamical (NEMO-OCE in a ORCA2 configuration [19]), a sea ice (NEMO-ICE) and a biogeochemical model (NEMO-PISCES) to simulate the ocean at global scale. The biogeochemical component, named Pelagic Interaction Scheme for Carbon and Ecosystem Studies (PISCES-v2), is an intermediate complexity marine ecosystem model simulating the main biogeochemical cycles and the first trophic levels. Four main living compartments are simulated (two phytoplankton and two zooplankton functional groups) as well as detritus (including two size classes, small and large organic particulate carbon (POC)), nutrients and O₂ [20]. As an extensive description and validation of the model has been published in [20], therefore only the relevant parts of the model are recalled here.

In the model, air–sea O₂ exchange is computed using the commonly used Wanninkhof's parameterization [21]. O₂ concentration in the ocean depends mostly on photosynthesis (O₂ production) and remineralization (O₂ consumption) with a O₂/C Redfield ratio equal to 1.34 [22]. Denitrification process (i.e. the use of the oxygen from nitrate (NO₃) when O₂ concentration lower than 6 µM [23]) is explicitly simulated in PISCES.

The relative complexity of the O₂ cycle in the NEMO-PISCES model results in good abilities to (i) represent volumes of low oxygenated waters at global scale [2] (figure C1(b)), and (ii) simulate a reasonable mean state compared to other global biogeochemical models [24, 25] used in CMIP. These good performances explain our choice of using the NEMO-PISCES model for the present study.

We used the PISCES model forced offline by an ocean dynamical seasonal climatology as in [20]. The model spatial horizontal resolution is 2^∘ by 2^∘ cosφ (where φ is the latitude) with a meridional resolution refined to 0.5^∘ near the equator. The vertical grid is composed of 31 layers of widths varying between 10 m near the surface to 500 m at depth. Prior to our experiments, we have spun up the biogeochemical model for 2000 years, starting from observed climatologies for nutrients [26] and from uniform low values for ecological tracers, until quasi-equilibrium was reached. After this spin-up, 100 additional years were simulated in two different configurations: (1) a control simulation (CTRL) without AO perturbation, (2) a simulation with an added AO (OXYBIO). A third configuration, named OXYNOBIO, simulates AO as a conservative passive tracer regardless of the oxygen cycle, with the exact same source of AO than the one used in OXYBIO. The aim of simulating the OXYNOBIO tracer is to track the transport of AO over a span of 100 years.

The experimental scenario used for the AO injection was based on the theoretical economical viability and industrial feasibility of offshore implantation of hydrogen Production Units (PUs). Offshore locations have been gauged considering socioeconomic and technological data mainly in line with the forecasts from the International Energy Agency [27]. These theoretical locations were provided by the LHYFE company [28] which is a producer and supplier of green hydrogen. This idealized scenario was designed to give an upper bound of an AO injection by assuming a global adoption of hydrogen as a key component of renewable energy supply chains by 2030, which corresponds to the highest level of plausible green hydrogen use.

For this study, the assumed number of hydrogen plants installed over the next 50 years follow roughly the same deployment than the offshore oil platforms between 1975 and 2020: approximately 20 000 platforms. A scenario allowing the deployment of 17 936 PUs has therefore been defined. The PUs have been distributed within 304 selected model grid points (each covering about 400 000 $\mathrm{km}^2$, called hereafter 'stations'), resulting in an average of 59 PUs per model grid point. These stations were distributed along the world's coasts (figure 1) according to the combination of two criteria: (i) an average wind speed preferably between 12–17 $\mathrm{m\,s}^{-1}$ (optimal operation of the wind turbines) or between the minimum and maximum wind speed needed to operate the wind turbine (5–25 $\mathrm{m\,s}^{-1}$ [27, 29]), and (ii) a population density that is presumed to exceed 2000 people $\mathrm{km}^{-2}$ (related to the size of the energy market), or display a constant or even positive natural evolution of the population by 2050 [27, 30]. As a result, the locations of the PUs are not necessarily positioned near known OMZs. Finally, we chose an injection depth of 150 m, except in areas with shallower bathymetry where the AO was injected at the greatest depth possible. Indeed, injection at a depth exceeding 150 m (e.g. close to the core of the OMZ) would have been more costly and would therefore not be economically viable. Nevertheless, the depth chosen here remains close to the upper limit of the OMZ [31] (figure C1).

Zoom In Zoom Out Reset image size

Currently, the green hydrogen production can reach or even exceed 129 000 t H₂ $\mathrm{y}^{-1}$, depending on the PU [28]. During the production of hydrogen by electrolysis, the H₂/O₂ ratio is 2:1 [32]. This amounts to 8 kg of O₂ per 1 kg of H₂ produced. Thus, according to LHYFE's forecasts, each of the PUs would produce 70 000 t O₂ $\mathrm{y}^{-1}$ (figure 1(a)), i.e. 39.2 $\mathrm{Tmol\,y}^{-1}$ globally, which corresponds to ${\backsim}$0.86% of the overall biogenic O₂ production in the CTRL simulation between 0–150 m (4.56 $\mathrm{Pmol\,y}^{-1}$).

It is important to note that the idealized setup of our experimental design leaves aside engineering injection details (discussed in section 3.5) that may impact the efficiency of the AO injection.

After 100 years of AO injection (OXYBIO run), the global O₂ inventory increased by 0.07% with respect to the CTRL run (estimated to be 233.1 Pmol, a value close to that obtained from the observations by [8]: 227.4 Pmol). Yet, the AO injection reduced the volume of OMZs by between 1.1% and 2.4%, depending on the O₂ concentration threshold used (table 1). This suggests that the overall evolution of the OMZ volume is weakly sensitive to the chosen threshold. However, the global average hides significant regional disparities.

Table 1. The proportion of the volume of each OMZ to the overall volume of OMZs (volume contribution) at the last simulated year of the CTRL run and averaged between 200 and 600 m depth, and change (OXYBIO minus CTRL) in volume of the oxygen-poor water (generically referred to as OMZs) after 100 years of AO injection. The volumes are defined by oxygen thresholds $80\,\mathrm{mmol\,m}^{-3}$, $50\,\mathrm{mmol\,m}^{-3}$ and $5\,\mathrm{mmol\,m}^{-3}$. Free fields (*) represent the absence of water masses with O₂ concentrations below $5\,\mathrm{mmol\,m}^{-3}$ and bold indicates values for OMZs at global scale (sum of OMZs).

OMZ	$80\,\mathrm{mmol\,m}^{-3}$		$50\,\mathrm{mmol\,m}^{-3}$		$5\,\mathrm{mmol\,m}^{-3}$
%	Change	Volume contribution	Change	Volume contribution	Change	Volume contribution
Global	−2,4	100	−1,1	100	−1,1	100
Atlantic Ocean	−13,3	6,6	−34,6	4	0	2,4
Pacific Ocean	−2,3	81,3	−0,6	86,1	−2	94,6
Bay of Bengal	+5,5	8,1	+11,6	7,9	+25,3	3
Arabian Sea	−3,4	4,1	−6,1	2	*	*

The OMZs that shrank the most are the least extended (figure 2(a), table 1): in the Atlantic Ocean and the Arabian Sea (AS), that shrank by 13.3%–35.6% and 3.4%–6.1% respectively ($80\,\mathrm{mmol\,m}^{-3}$ and $50\,\mathrm{mmol\,m}^{-3}$ thresholds). The smallest shrinkage (0.6%–2.3%) is that of the Pacific Ocean OMZ (table 1). The most striking consequence of the AO injection is simulated in the Bay of Bengal (BoB) where the OMZ expands counterintuitively, its volume increases by 5.5%–11.6% for the $80\,\mathrm{mmol\,m}^{-3}$ and $50\,\mathrm{mmol\,m}^{-3}$ thresholds, and even 25.3% for the $5\,\mathrm{mmol\,m}^{-3}$ threshold (table 1 and figure 2(b)).

Zoom In Zoom Out Reset image size

Contrary to the global assessment, the volume variation of regional OMZs depends strongly on the threshold that is considered to define them. The harshest thresholds specific to each region allow us to determine the evolution of the OMZ core. Thus, excluding the Pacific OMZ, the volume change is the most relevant in the OMZ core (i.e. a larger volume variation for the threshold of $50\,\mathrm{mmol\,m}^{-3}$ in the Atlantic Ocean and $5\,\mathrm{mmol\,m}^{-3}$ in the AS and BoB than for $80\,\mathrm{mmol\,m}^{-3}$, table 1). However, the evaluation of the structure of an OMZ requires to also consider its upper and lower boundaries (oxyclines) [31]. OXYBIO configuration simulated a slight deepening of the upper oxycline of all OMZs (figures 2(c)–(f)), more marked in the Atlantic Ocean. Thus, their thickness thins by an average of 42 m (for the Atlantic Ocean) and 13 m (for the AS) between 0 and 600 m while that of the BoB thickens by 11 m despite the substantial increase in O₂ between 100 and 200 m, which can exceed 100% (figure 2(e)). The shift of the horizontal boundaries position of the OMZs in addition to the deepening of the vertical boundaries, characterize the strongest shrinkage of the OMZs, namely in the Atlantic Ocean, and to a lesser extent in the AS (figures 2(c) and (f)).

A wide range of ecosystem processes interacting with the three-dimensional circulation determine the distribution of O₂ concentrations in the ocean [33]. The analyses in section 3.1 suggest a regional disparity in the mechanisms causing O₂ variation. $O_{2}^\mathrm{OXYNOBIO}$ (AO in the OXYNOBIO configuration) tracks exclusively the transport of AO without any air–sea exchanges and biogeochemical reactions. The correspondence (R) between the variation of $O_{2}^\mathrm{OXYNOBIO}$ and that of $O_{2}^\mathrm{(OXYBIO-CTRL)}$ (i.e. O₂ change) addresses the relative importance of biological activity and transport to O₂ variation in OXYBIO (figure 1 (b), equation (1) and appendix A).

$$\begin{eqnarray} &&R = \frac{\mathrm{Variance}\left(O_2^\mathrm{(OXYBIO-CTRL)}-O_2^\mathrm{OXYNOBIO}\right)}{\mathrm{Variance}\left(O_2^\mathrm{OXYNOBIO}\right)}. \end{eqnarray} \tag{ 1 }$$

R contrasts O₂ changes due to biological response (i.e. $O_{2}^\mathrm{(OXYBIO-CTRL)}$ minus $O_{2}^\mathrm{OXYNOBIO}$) to AO transport ($O_{2}^\mathrm{OXYNOBIO}$). Thus, in (figure 1(b)) the mechanisms are categorized as follow:

• (i)  
  R$\backsim$0 indicates the dominance of AO dynamic transport over biological activity as a mechanism for the variation of O₂,
• (ii)  
  R$\backsim$1 indicates that transport and biogeochemical processes are of similar importance to the O₂ changes. This suggests either a biological change concomitant with AO transport or transport of what result from remote biological response,
• (iii)  
  R$\gt$1 indicates the dominance of biological activity over dynamic transport to explain the variation of O₂ consecutive to AO injection.

Along the coastlines with injection sites, in the Atlantic Ocean and across the North Pacific Equatorial Current, the correspondence is high (R$\backsim$0). In the Atlantic Ocean, the OMZ shrank even where it is devoid of AO stations (figure 1(b)). Advective transport is therefore the driving factor feeding AO into this OMZ. Since this OMZ is included in the North Atlantic Central Water (figure C2) and has its top located in the Subtropical Mode Water (figure 2(c)), it could benefit from an efficient supply of AO. This OMZ would therefore be fed by the small cell-dominated by meridional surface transport as well as by the upper limb of the Atlantic Meridional Overturning Circulation (AMOC, figure C2). An exact quantification of this would require a dedicated study. In the South Pacific Ocean and the Costa Rica Dome, a combined effect of biogeochemical feedbacks and dynamic transport ($R\backsim$1) seems to explain the response to AO injection (figure 1(b)).

In the Indian Ocean, biogeochemical feedbacks dominate, especially in the BoB ($R\gt$1, figure 1(b)) where they appear to explain the counter-intuitive reduction of O₂ concentration in the OXYBIO simulation with respect to CTRL (figure 2(b)) as treated in section 3.3. It is noteworthy that the different biogeochemical responses between the AS and BoB is likely due to different levels of oxygen deficiency [34].

In the ocean, O₂ is consumed by biological respiration, including OM remineralization. When the waters become suboxic, oxygen from NO₃ is used for remineralization, a process called denitrification. OM remineralization is thus an internal source of 'recycled' fixed nitrogen (N) except in suboxic waters where denitrification is a sink of bioavailable nitrogen. In this section we investigate the potential cascading effects of the AO on the denitrification, remineralization, and Net Primary Production (NPP) (figures 3 and 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Denitrification is absent in the Atlantic Ocean OMZ (figure C3) due to relatively greater natural ventilation than in the other OMZs. The decrease (up to 10% in the Pacific Ocean) or even total inhibition (in the BoB) of denitrification in response to AO injection (figure 3(a)) reduces the NO₃ sink in initially (before AO injection) suboxic areas (figure 3, orange line). Note that at the surface, nutrients decrease except NO₃ which increases (figure 3(c)), while at the subsurface all nutrients increase (figures C4, C5 and 3(d)). This suggests that the reduction of denitrification is likely to increase the availability of NO₃ inducing the stimulation of NPP (figure 4(b)), thus increasing the exported production (figure 5(a)). As remineralization is constrained by OM [35–37], the additional supply of OM can thus amplify and maintain high remineralization leading to an increase in deep nutrient stocks [1]. Indeed, the seasonal thermocline of BoB is a major nutrient reservoir through the intense remineralization of exported OM [38]. However, by increasing NPP, nutrient uptake is amplified. Thus, the originally non-limiting nutrients may decrease until depletion [39] as in the BoB where NO₃ limitation shifts to phosphate limitation for nanophytoplankton and changes from silicate to phosphate limitation for diatoms (figure C6).

Zoom In Zoom Out Reset image size

In conclusion, the expansion of OMZ indicates an imbalance between oxygen supply and biological O₂ consumption, as in the case of BoB. Despite continuous oxygenated water injection and O₂ transport through ocean dynamics, biological activity can still drive an expansion of an OMZ in some regions. The suppression of the nitrate sink (by inhibiting denitrification) and enhancing nutrient recycling (through remineralization) particularly at the seasonal thermocline is the cause. Moreover, the decline of the denitrification results in a decrease in nitrogen fixation (N₂-fixation, roughly 50% in the vicinity of BoB) (figure C7). Landolfi's study [40] emphasized the vicious cycle induced by the spatial coupling of denitrification and N₂-fixation. Due to stoichiometric imbalance, The NO₃ loss through denitrification versus the nitrogen gain through the OM remineralization derived from N₂-fixation leads to a net loss of fixed nitrogen that further stimulates N₂-fixation. Counter-intuitively, the reduction in N₂-fixation in the OXYBIO simulation, through the ihnibition of denitrification, could thus contribute to an increase in fixed nitrogen.

The other consequence particularly notable is the nonlocal environmental effects of AO injection. The AO stimulates biological activities beyond the zones of O₂ change, especially in the Indian and tropical Pacific Oceans, except in the Northern Hemisphere where O₂ changes and biological activities tally (figure 1(b)). While denitrification changes are mainly confined to the continental margins (figure 3(a)), remineralization changes are remotely triggered (figure 3(b)). Whether it is related to aerobic respiration (figure 3(b)) or NPP (figure 4(b)), biogeochemical changes follow the pathways of the advective currents (i.e. circulation of the subtropical gyres) (figures 1(b) and 4). In the Eastern Pacific Equatorial zone, phytoplankton is mainly Fe-limited in the CTRL simulation (figure C6). Thus, the gain in NO₃ (figures 3(c) and (d)), due to reduced denitrification in this zone, would feed the initially NO₃-limited phytoplankton communities in the areas surrounding the Pacific cold tongue and along the Central American coasts (figures C6(a) and (c)), through lateral surface transport. As a result, NPP increases north of the equator (figure 4(b)) along the Central American coasts (from the Gulf of Panama to the Mexican coast), where there are no injection stations (figure 1). This local NPP increase would trigger O₂-consuming remineralization (figure 3(b)). Thus, the Costa Rica Dome, an eddy structure particularly influenced by the Costa Rica Coastal Current from the south, would be subject to exacerbated denitrification (up to 5%, figure 3(a)) due to this stimulation of remineralization bordering the dome. In the Tropical South Pacific Ocean, the edges of the denitrification zones (where it is lowest (figure C3), indicated by the orange line in figure 3(a)), show the largest decrease in denitrification, by 4%. Beyond these edges, the NPP remote response (rising to +3% off Chile (figure 4(b)) is triggered by the new supply of NO₃ in zones that can be limited in NO₃ (figure C6), which could explain the lack of NPP response near coasts where phytoplankton is mainly limited by Fe and Si (figure C6). Hence, these biogeochemical changes would be transported along the South Equatorial Current, feeding the subtropical gyre through the lateral transport [41].

Another striking region is the Indian Subtropical oligotrophic gyre (figure 4(a)), where the biological response is opposite to all other basins (figures 3(b) and 4). In BoB, the strong increase in NPP results in increased nutrient consumption (figure C4). This drop in nutrients could lead to a reduction in the supply of nutrients to the Indian Subtropical gyre through the South Java Current and explain the NPP decrease. In such a 'regenerative loop' regime [42] where most production is driven by recycled nutrients [43], the small reduction in OM remineralization may also explain the decrease in NPP in a chain-reaction.

Fertilization, whether natural [44] or artificial [45], can favor exported production through the development of ecosystems prone to exceptional blooms. These ecosystems are based on a biomass dominated by large cells (e.g. diatom [46–48]) and a NPP provided jointly by small and large cells [49]. The AO injection-induced denitrification inhibition does not appear to lead a change in the phytoplankton community despite NO₃ enrichment (figure C8). Besides phytoplankton community composition, estimation of NPP is important for inferring the additional sequestered carbon. AO injection-induced NPP increase (+0.7%) is close to the magnitude of the decrease predicted by the low-emission scenario SSP1-2.6 ($-0.56\pm4.12$%) but is less than the decrease predicted by the high emission scenario SSP5-8.5 ($-2.99 \pm 9.11$%) by 2100 [11], which assumes a potentially negligible effect on the overall air–sea CO₂ flux.

Because of the ability of artificial oxygenation to increase nitrogen concentration, it can also be suggested as an indirect means of fertilization by industry to stimulate the carbon pump. Yet, iron fertilization experiments revealed inconsistent export of production despite increased NPP [50]. In our study, despite the joint increase in NPP and exported production (figures 4(b) and 5(a)), the e-ratio decreases between 1%–14% depending on the region (figure 5(b)), which indicates a decrease in the efficiency of carbon export. Thus, due to a regionally variable duality between remineralization processes and OM export beyond injection sites, the effect of AO on carbon storage should not be overlooked in future studies.

Several sources of uncertainty are possible when analyzing the effects of AO injection on the global ocean. The deployment scenario used in this study may represent one of them and may affect the results. For example, the choice of injection depth is constrained by structural robustness and cost. However, at 1000 m the oxygen utilization rate is 10 to 50 times lower than at 100 m [10, 51–54]. Hence, depth is an important causal parameter in the biogeochemical and ecosystem response to AO injection. The coarse resolution of the model used (2^∘ resolution) prevents us from documenting technologies required to inject AO, which come under sub-grid processes but may affect AO efficiency [1, 15]. Yet, we can infer that AO should be undersaturated at the injection sites (figure C9). Thus, O₂ saturation at injection sites (with no technical limitation taken into account) would allow the full dissolution of the AO. In addition, the spatial resolution implies process-related biases below the mesoscale spectrum [55]. Regional modeling of the AS highlighted the importance of eddy-induced advection on vertical and lateral O₂ transport [56]. Increasing the resolution improves the explicit representation of eddy mixing that would be involved in the biased volume representation of OMZs in noneddy resolving ocean models. Thus, eddy mixing may be crucial in shaping marine habitats of oxygen-sensitive species [57].

Numerical modeling of OMZs and their ecosystems is an important socioeconomic issue because they are important fishing areas [58–60]. However, the current ocean microbiome modeling remains simple compared to the complexity that emerges from observations [61]. Also, AO injection showed that the evolution of carbon transfer to the deep ocean [62] as well as the fate of the OMZs, and consequently of the fisheries [58–60], require the analysis of phytoplankton and zooplankton communities. In our study, the microbial recycling mechanisms are parameterized lacking an explicit representation of bacterial communities despite the influence of O₂ on controlling the availability of fixed nitrogen in the ocean [63, 64]. Regionally, explicit modeling of some ecosystems can improve the simulation of these processes, such as the diazotrophic community in the Tropical Pacific [65], especially since the anaerobic to aerobic transition in the ocean can cause strong negative feedback on nitrogen fixation [64]. Ultimately, ecosystem diversity differs according to the model used, and uncertainties in their applications remain [3, 61, 66, 67] and limit management or geoengineering decisions that are based solely on the use of models [68].

Despite these uncertainties, our results describe changes in biogeochemical processes in line with published regional studies on the effects of changing oxygen levels (e.g. in the Arabian Sea [57]). In the Pacific ocean, oxygenation through the modification of modeled vertical water exchange have also ultimately resulted in a deoxygenation and remote changes in biological activity [1]. However, given that biogeochemical processes tend towards a new equilibrium state after 20 to 40 years of simulation, as in BoB (not shown), we do not prejudge the effects of extending the experiment, nor of stopping the injection, either during or after the 100 year experiment, which could potentially lead to a global and regional new state.

Vinogradova NT, Ponte RM, Stammer D. 2007 Relation between sea level and bottom pressure and the vertical dependence of oceanic variability. Geophys Res Lett. 34. https://doi.org/10.1029/2006GL028588. Bingham RJ, Hughes CW. 2008 The relationship between sea-level and bottom pressure variability in an eddy permitting ocean model. Geophys Res Lett. 35. https://doi.org/10.1029/2007GL032662. Talley LD, Pickard GL, Emery WJ, Swift JH. 2011 Descriptive Physical Oceanography. Ed Elsevier. https://doi.org/10.1016/C2009-0-24322-4. Pollard, R. T., & Pu, S. 1985. Structure and circulation of the Upper Atlantic Ocean northeast of the Azores. Progress in Oceanography, 14(C), 443–462. https://doi.org/10.1016/0079-6611(85)90022-9.