Controlling factors of the climatological annual cycle of the surface mixed layer oxygen content: A global view

Wang, Zhankun; Garcia, Hernan E.; Boyer, Tim P.; Reagan, James; Cebrian, Just

doi:10.3389/fmars.2022.1001095

ORIGINAL RESEARCH article

Front. Mar. Sci., 27 October 2022

Sec. Marine Biogeochemistry

Volume 9 - 2022 | https://doi.org/10.3389/fmars.2022.1001095

Controlling factors of the climatological annual cycle of the surface mixed layer oxygen content: A global view

Zhankun Wang^1,2*

Hernan E. Garcia³

Tim P. Boyer³

James Reagan³

Just Cebrian^1,2

¹National Oceanic and Atmospheric Administration’s National Centers for Environmental Information, Oceanographic and Geophysical Science and Services Division, Stennis Space Center, MS, United States
²Northern Gulf Institute, Mississippi State University, Stennis Space Center, MS, United States
³National Oceanic and Atmospheric Administration’s National Centers for Environmental Information, Oceanographic and Geophysical Science and Services Division, Silver Spring, MD, United States

The annual cycle of global dissolved oxygen content (O₂C) and mean oxygen concentration in the surface mixed layer are estimated using monthly climatological oxygen fields from the World Ocean Atlas 2018 (WOA18). The largest seasonal variability in the mixed layer O₂C occurs in the extra-tropics between 30° and 70° latitude of each hemisphere. A global view of the role of entrainment, air-sea flux, and biological activity in controlling the oxygen content/concentration annual cycle in the mixed layer is determined using an oxygen mass balance model. Based on the relative percentage from the mass balance model, entrainment is only a significant driver (contributing to 20-40% of the total changes) from mid-fall to early spring when the mixed layer deepens and transfers oxygen to deeper waters. Both the air-sea oxygen flux and biological activity show strong annual cycles and play critical roles in the annual cycle of O₂C in the mixed layer. Air-sea oxygen flux is ingassing from late fall to early spring and outgassing between late spring and early fall. It is a substantial factor throughout the year and controls 40-60% of the oxygen changes in most months. Biological activity is a net source (production) in the spring and summer and a net sink (consumption) in the late fall and winter for the mixed layer oxygen content. Biological activity plays a more important role during spring/summer (40-60%) than that during fall/winter (10-30%) in controlling the overall oxygen change in each month. The model estimates a mean value (±SD) of 3.06±1.61 mol C m^-2 yr^-1 and a total of 863.7±73.8 Tmol C yr^-1 for the global annual net ocean community production (ANCP) between 60°S and 60°N latitude, which are in fairly good agreement with previous studies.

1 Introduction

Understanding the controlling factors of the seasonal patterns of dissolved oxygen in the upper global ocean is important for oxygen dynamics and deoxygenation research in the upper global ocean (Garcia et al., 2005a; Garcia et al., 2005b) because the seasonal patterns typically need to be removed to focus on the long-term deoxygenation trends (Ito et al., 2017). Seasonal changes in the surface mixed layer are essential in understanding the dynamical ocean and climate, such as the air-sea oxygen fluxes, ocean stratification, mixing and biological activity, because the surface mixed layer constantly exchanges energy and gas with the atmosphere. In the surface mixed layer, oxygen is supplied through air-sea gas exchange (Najjar and Keeling, 2000; Garcia and Keeling, 2001; Bushinsky et al., 2017) and from photosynthesizing marine plants (Oschlies et al., 2018). The oxygen is then redistributed via entrainment, circulation, diffusion/mixing, consumption or respiration by marine life, and by chemical processes such as reduction-oxidation.

Air-sea exchange is likely the leading physical process determining oxygen annual cycle in the mixed layer and has been extensively studied (e.g. Najjar and Keeling, 2000; Garcia and Keeling, 2001; Emerson and Bushinsky, 2016; Bushinsky et al., 2017). Air-sea flux is conducted through diffusive oxygen transfer (Najjar and Keeling, 2000; Wanninkhof, 2014) and bubble injection (Liang et al., 2013; Emerson and Bushinsky, 2016). The diffusive transfer can be estimated as a function of wind speed and the difference between the measured O₂ and the oxygen saturation concentration (Najjar and Keeling, 2000; Liang et al., 2013). At the base of the mixed layer, the O₂ in the mixed layer mixes with deeper layers through entrainment (deepening of the mixed layer), turbulent mixing induced vertical diffusion (Ito et al., 2022), and upwelling/downwelling. Horizontally, O₂ is redistributed by diffusion and geostrophic and Ekman advection (e.g. Bushinsky and Emerson, 2018). Large scale circulations and eddies can also cause exchanges and mixing of water masses with different oxygen concentrations. Furthermore, the variability in O₂ content can be influenced by other physical processes. For example, Ito et al. (2019) have found that the dominant mode of O₂ variability in the North Pacific is significantly correlated with the Pacific Decadal Oscillation (PDO), which is driven by two primary mechanisms: vertical movement of isopycnals determining the oxygen variations in the tropics, and changes in subduction driving the variability in the subtropics. They also indicate that climate-change driven biological processes may reinforce the physically driven O₂ change in the tropics. Tides and atmospheric forcing can also impact O₂ content in the intermediate waters of the subarctic North Pacific (Andreev and Baturina, 2006).

Biochemical processes are typically difficult to measure directly compared to physical processes in O₂ annual cycle studies. One approach for determining biogeochemical processes is to use a mass balance model by studying the oxygen annual cycle in the mixed layer and calculating the residual oxygen flux by subtracting air-sea flux, entrainment, diffusion and transport terms from the oxygen content change (e.g. Bushinsky et al., 2017; Yang et al., 2017; Yang et al., 2018). When dominant physically-induced oxygen fluxes are eliminated, the residual is a rough estimate of the net “biochemical” oxygen production. Using a stoichiometric relationship between oxygen evolved and net community production (NCP), oxygen budget can then be used to determine the annual cycle in NCP (e.g. Keeling and Shertz, 1992; Jin et al., 2007; Bushinsky et al., 2017; Yang et al., 2018). Numerical models have also been used to quantify the roles of biochemical and physical processes on seasonal oxygen variations (e.g. Jin et al., 2007; Manizza et al., 2012).

In this study, we first analyze the seasonal variations of the O₂C and mean oxygen concentration in the surface mixed layer for the global ocean using objectively-analyzed monthly climatology fields from the World Ocean Atlas 2018 (WOA18, Boyer et al., 2018; Locarnini et al., 2018; Garcia et al., 2018). We then, from a global view, examine the role of three dominant processes controlling the O₂C annual cycle in the mixed layer. A mass balance model is utilized to study the role of entrainment, air-sea flux and biological activity in the oxygen content seasonal change. Similar oxygen mass balance models have been used for oxygen dynamics studies at fixed stations (e.g. Emerson, 2014) and for regional oceans (e.g. Bushinsky et al., 2017; Yang et al., 2018). Utilizing the global oxygen monthly climatology fields and the newly-developed global mixed layer depth fields from the WOA18, we apply the mass balance model to the global ocean in this study. A better characterization of the nature, variability and drivers of the annual cycle of the global oxygen content in the mixed layer will improve our understanding and predictability of interannual variability and long-term trends (e.g. deoxygenation) in oxygen content in the global ocean (Garcia et al., 2005b; Jin et al., 2007).

2 Data and methods

2.1 Climatological oxygen fields and WOA18

Oxygen fields are from the WOA18 objectively-analyzed fields on one-degree longitude/latitude grids (Locarnini et al., 2018; Garcia et al., 2018). The WOA18 climatologies are “all-data” climatologies that are created using all available quality-controlled data regardless of year of observation collected from 1960 to 2017 on 102 standard depth levels from the surface to 5,500 m (See Table 3 in the WOA document for more details about the standard levels, https://www.ncei.noaa.gov/sites/default/files/2022-06/woa18documentation.pdf) . WOA18 contains annual, seasonal, and monthly climatologies and related statistical fields: the number of observations, standard deviations, and standard errors for each grid and depth. Oxygen data used to create the climatologies are from the discrete water sample dataset in the World Ocean Database (WOD18, Boyer et al., 2018), obtained by chemical Winkler titration methods collected between 1960 and 2017. All data are analyzed in a consistent manner by a series of quality control (QC) procedures. The QC includes duplicate elimination, range and gradient checks, statistical checks and subjective flagging. The objective analysis/mapping uses a three-iteration scheme of the type described by Barnes (1964) with a Gaussian weight function and variable influence radii of 892, 669, and 446 km (Locarnini et al., 2018). More details on data sources and distribution, data QC and objective analysis can be found in the WOA documents (Garcia et al., 2018; Locarnini et al., 2018; Zweng et al., 2018). We follow the WOA season definition in our analysis (namely, spring: March-May; summer: June-August; fall: September-November; and winter: December-February for the Northern Hemisphere (NH)).

2.2 Climatological mixed layer depth

A new 1°×1° spatial latitude-longitude resolution global climatology of the mixed layer depth (MLD) based on individual profiles from WOD18 is constructed as part of WOA18. For each cast in WOD18 that has concurrent measurements of temperature and salinity, potential density (with a reference level of 0 m) is calculated using the Gibbs-SeaWater Oceanographic Toolbox (McDougall and Barker, 2011), which contains the Thermodynamic Equation of State - 2010 (TEOS10; IOC, SCOR and IAPSO, 2010) subroutines. While many different criteria exist for estimating the MLD, we choose to follow the density threshold criterion (Δσ_θ=0.125kg m^-3) which estimates MLD by approximating the depth for which potential density changes 0.125kg m^-3 from a near-surface value at 10 m depth. We choose this density criteria because it reduces underestimations when compared to smaller thresholds (e.g., 0.03 kg m^-3, Toyoda et al., 2017) and is frequently used in climate-related research where MLD is often averaged over monthly and longer time periods (Kara et al., 2000; Thomson and Fine, 2003). Our procedures follow a similar method described in de Boyer Montégut et al. (2004) which estimated MLD from individual profiles with most data at observed levels. The reference depth is selected at 10m to minimize the strong diurnal cycle in the top few meters of the ocean (de Boyer Montégut et al., 2004). Therefore, our shallowest mixed layer depth is 10 m.

2.3 Oxygen content calculation

WOA18 monthly O₂ fields on a 1°×1° grid are used to calculate the monthly O₂C in the surface ML for each 1°×1° grid of the global ocean,

\begin{array}{l} O 2 C = A \int_{0}^{h} O_{2} d z = A h O_{2 m} & (1) \end{array}

where A is the area in m² of each 1°×1° latitude-longitude grid box, dz is half the distance between the next shallowest level and the current level plus half the distance between the next deepest level and the current level in the WOA (i.e., thickness of each depth layer), and h is the MLD from section 2.2. $O_{2 m} = \frac{1}{h} \int_{0}^{h} O_{2} d z$ is the average oxygen concentration in the surface mixed layer.

We define zonal and basin oxygen content as the sum of the O₂C integrated along a 1° latitude belt over the global or basin scales. To compare values among basins with different areas, average zonal and basin O₂C are used, which are calculated from the integrals divided by the areas of the study regions.

2.4 Oxygen mass balance model

Seasonal changes in the mixed layer O₂ concentration result mainly from air-sea oxygen exchange, entrainment, and biochemical activity within the mixed layer. To gain a better understanding of the role of these three controlling factors, we use an oxygen mass balance model following a schematic like Emerson et al. (2008); Emerson and Stump (2010); Bushinsky and Emerson (2015) and Yang et al. (2017); Yang et al. (2018). In this model, O_2m variability in the ML is balanced with air-sea interface O₂ exchange flux F_A-W, entrainment F_E and biological activity J_O2. Horizontal advection/diffusion and vertical diffusion are assumed to be small as discussed in Jin et al. (2007); Bushinsky and Emerson (2015) and Yang et al. (2017). We use oxygen produced by “biological activity” to refer to J_O2 in our following discussions, but keep in mind that it might also contain other terms, such as changes due to chemical processes, transport, diffusion, and errors.

The equation for mean oxygen concentration changes in the surface mixed layer may be written as (Emerson et al., 2008; Emerson and Stump, 2010)

\begin{array}{l} \frac{h \partial O_{2 m}}{\partial t} = F_{A - w} + F_{E} + J_{O 2} & (2) \end{array}

where ∂ indicates a partial derivative. A detailed explanation of how F_E and F_A-W are calculated can be found in Section 2.4.1 and 2.4.2 below respectively and more details in Emerson and Stump (2010); Bushinsky and Emerson (2015), and Yang et al. (2017). The term on the left-hand side of Eq. (2) represents the O₂C change related to the O_2m concentration change (see Supplemental Material Section 1). The “biological” oxygen production, J_O2, is calculated as the residual $J_{O 2} = \frac{h \partial O_{2 m}}{\partial t} - F_{A - w} - F_{E}$ .

2.4.1 Entrainment F_E

The entrainment of oxygen can be estimated as (Emerson et al., 2008; Ren and Riser, 2009)

\begin{array}{l} F_{E} = - w_{e} δ O_{2} = - w_{e} (O_{2 m} - O_{2 b}) & (3) \end{array}

Where w_e is the entrainment velocity, and δO₂ is the oxygen difference between the average O₂ concentration within the mixed layer, O_2m, and the O₂ of the first layer below the base of the mixed layer, O_2b, in the WOA. The entrainment velocity is estimated as (Stevenson and Niiler, 1983; Ren and Riser, 2009),

\begin{array}{l} w_{e} = H (\frac{\partial h}{\partial t} + \nabla \cdot h \vec{v}) & (4) \end{array}

Where h is the mixed layer depth, t is for time and $\vec{v}$ is the horizontal velocity including Ekman velocity and geostrophic velocity. We calculate the monthly climatological horizontal velocity from the Simple Ocean Data Assimilation (SODA) reanalysis dataset (https://dsrs.atmos.umd.edu/DATA/soda3.4.2/REGRIDED/) by averaging all the data between 1980 and 2019 for each month. Upwelling/downwelling is considered in the entrainment velocity. H is the Heaviside unit function

\begin{array}{l} H (x) = {\begin{cases} 1, x \geq 0 \\ 0, x < 0 \end{cases} & (5) \end{array}

Only entrainment (w_e ≥ 0) is considered and detrainment (w_e< 0) is set to zero because the water flowing out from the mixed layer has the same O₂ concentration as the water in the mixed layer and will not affect the mixed layer O₂ concentration.

2.4.2 Air-sea Oxygen Flux F_A-W

The total air-sea oxygen flux F_A-W can be estimated by the sum of the diffusive oxygen flux F_S and bubble induced air-sea flux (Liang et al., 2013; Emerson and Bushinsky, 2016; Bushinsky et al., 2017),

\begin{array}{l} F_{A - w} = F_{s} + β (F_{C} + F_{p}) & (6) \end{array}

where F_C is the oxygen flux due to small bubbles that completely collapse, F_p is the flux due to large bubbles that partially collapse. Bubble processes can be estimated using different bubble models (e.g. Woolf and Thorpe, 1991; Woolf, 1997; Stanley et al., 2009; Liang et al., 2013). Through comparison of noble gas and N₂ measurements, Emerson and Bushinsky (2016) suggest that the bubble model of Liang et al. (2013) is preferred for air-sea O₂ flux calculation because it includes mechanisms for both large and small bubbles. We will use the same bubble model of Liang et al. (2013) in the present study. β=0.37 is a correction factor for the bubble-mediated oxygen flux (Emerson et al., 2019). In all the air-sea flux components, positive is to the ocean (ingassing) and negative is to the atmosphere (outgassing).

F_s can be calculated following the formula in Bushinsky et al. (2017) and Liang et al. (2013) and considering the adjustments for sea level pressure and skin temperature, respectively (Najjar and Keeling, 1997; Najjar and Keeling, 2000)

\begin{array}{l} F_{S} = - 1.3 \times 10^{- 4} U_{a}^{*} {(\frac{S c_{O_{2}}}{660})}^{- 0.5} (O_{2} - O_{2}^{s a t} + δ_{p} + δ_{T}) & (7) \end{array}

Where $U_{a}^{*} = \sqrt{C_{d} U_{10}}$ is the air-side friction velocity (Large and Pond, 1981; Emerson and Bushinsky, 2016). C_d is the drag coefficient, which is a function of the 10m wind speed U₁₀ and assuming the stability effect can be ignored (Kara et al., 2007). Sc₀₂ is the Schmidt number (Najjar and Keeling, 2000; Wanninkhof, 2014). $O_{2}^{s a t}$ is the oxygen saturation (Garcia and Gordon, 1992). The 10m wind speed data is averaged from the Cross-Calibrated Multi-Platform (CCMP) 6 hourly high-resolution data between 1988 and 2020 (http://www.remss.com/measurements/ccmp/). CCMP is a level 3 blended wind product which uses satellite, in situ data, and numerical modelling simulations (Wentz et al., 2015). Both O₂ and $O_{2}^{s a t}$ are from the WOA18. δ_p is the adjustment for sea level pressure variation and δ_T is the adjustment for the skin temperature (Najjar and Keeling, 2000). Both sea level pressure and skin temperature monthly climatologies were calculated from the European Centre for Medium-range Weather Forecasts (ECMWF) Reanalysis 5^th Generation data (ERA5) (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels-monthly-means?tab=form).

F_C may be calculated as (Liang et al., 2013; Emerson and Bushinsky, 2016),

\begin{array}{l} F_{C} = 5.56 \times {(U_{w}^{*})}^{3.86} X_{O 2} & (8) \end{array}

Where $U_{w}^{*} = \sqrt{P_{a} / P_{w}}, U_{a}^{*}$ is the ocean-side friction velocity and P_a is the air density and p_w is the water density. X_O₂ is the mole fraction of oxygen in the atmosphere.

F_p may be estimated using (Liang et al., 2013; Yang et al., 2017; Emerson et al., 2019; Huang et al., 2022):

\begin{array}{l} F_{p} = - 5.5 {(U_{w}^{*})}^{2.76} {(\frac{S c o_{2}}{660})}^{- 2 / 3} [O_{2} - (1 + Δ P) O_{2}^{s a t} + δ_{p} + δ_{T}] & (9) \end{array}

Where $Δ P = 1.52 {(U_{w}^{*})}^{1.06}$ is the increase in oxygen saturation within large bubbles due to the wind induced deepening effect (Liang et al., 2013). The units of the symbols used in the equations of section 2.4.2 are shown in Table 1.

TABLE 1

Table 1 The meaning of symbols used in Section 2.4.2 and their units.

3 Results and discussion

3.1 Zonal and seasonal variation of MLD

Most seasonal variability of the MLD occurs in the extra-tropics and subpolar regions within the 30° to 70° latitude belt of the Northern Hemisphere (NH) and Southern Hemisphere (SH) (Figure 1). In the NH, the largest (deepest) MLD is found in the Atlantic Ocean (Figures 1B, E) with depths greater than 1000m, which occurs around 60°N in March. The mixed layer starts to become shallower in April and rapidly decreases to less than 40 m on average in June. The average MLD remains shallower than 50 m between June and September for the global ocean and in each of the major ocean basins. The deepening of the ML starts in September and continues until it reaches the maximum value once again in March. The MLD in the SH shows a similar annual cycle, but differs by six months. The largest annual variation in the SH occurs in the Southern Ocean. The basin MLD climatology (Figures 1B–D, G–I) shows basin scale differences among the three basins. A full-cycle deep ML belt is formed in the Southern Ocean during SH winter between 40 and 60°S (spatial plots, not shown). The belt is shallowest in the Atlantic Ocean and deepest in the Pacific Ocean. The center location of the Southern Ocean ML belt is around 40° S in the Indian Ocean (Figure 1C) and shifts to 51°S in the Pacific Ocean (Figure 1D). The belt is much shallower in the Atlantic Ocean centered around 53°S. The monthly time series of average MLD climatologies between 30° and 70° in Figures 1F–I also shows a strong annual cycle in the subtropics and subpolar regions.

FIGURE 1

Figure 1 Average climatology of mixed layer depth (MLD) in meters. Zonally averaged monthly MLD for the (A) World, (B) Atlantic, (C) Indian, and (D) Pacific Ocean. Same colorbar is used for panels (A–D) as shown below them. The contour lines are 20, 30, 40, 50, 75, 100, 150, 200, 300, 400, 500, 600 and 700m. (E) Zonal and annual average of MLD as a function of latitude for the World Ocean and the three major ocean basins. Zonally averaged monthly time series estimate of the mixed layer depth for the subtropics and mid-latitudes between 30° and 70° for both Northern Hemisphere (NH, solid black) and Southern Hemisphere (SH, dashed black) for (F) World, (G) Atlantic, (H) Indian, and (I) Pacific Ocean, respectively. The shadow areas in (F–I) represent the data variability between one standard deviation (SD) above and one standard deviation below the average. The standard deviation is estimated based on the zonal averages shown in the plots (A–D) above.

3.2 Zonal and seasonal variation of O₂C and O_2m in the ML

Similar to the MLD seasonal variations, the largest seasonal changes in O₂C also occur in the extra-tropics and subpolar regions between 30° and 70° of the NH and SH (Figure 2). In the NH, the largest seasonal peak (>110.0 mol/m²) occurs between February and April and centered near 60°N (Figure 2A), which is when and where the MLD is deepest in the NH (Figure 1). The smallest O₂C (<10.0 mol/m²) occurs during summer between July and September in the NH. In the SH, the largest seasonal peak occurs approximately in September and centered at 55°S in the southern Pacific Ocean. The magnitude of the SH peak (~90.0 mol/m²) for the global ocean is slightly smaller than that of the NH. The center locations of the Southern Ocean O₂C peaks in the three basins (Figures 2B, C) are consistent with the center locations of the deepest MLD shown in Figure 1. In the tropical and subtropical regions, the seasonal variations are smaller compared to those in high latitudes. As would be expected, the seasonal and spatial variations of the ML O₂C for both global and the three major basins (Figure 2) are highly correlated with the seasonal and spatial distributions of the MLD (Figure 1) by comparing the two figures. In the northern Atlantic and southern Pacific, the mixed layer has larger seasonal variations and can hold 10+ times more O₂ content in the winter than that in the summer. Although MLD plays a dominant role in determining the peak O₂C in the winter, the O₂ concentration might be also larger in the winter than that in the fall due to enhanced air-sea flux and thermal-induced solubility increase (Najjar and Keeling, 1997; Najjar and Keeling, 2000; Garcia et al., 2005a). The winter O₂C peak in the extra-tropics and subpolar results from the combination of mixed layer deepening and O₂ solubility increase. The O₂C seasonal variation in the northern Pacific is much smaller than that in the northern Atlantic. The annual averaged O₂C in Figure 2E also shows huge spatial differences among different basins and the global ocean.

FIGURE 2

Figure 2 Climatology of area-averaged O₂C in the surface mixed layer (mol/m²). Zonally averaged monthly O₂C in the mixed layer for the (A) World, (B) Atlantic, (C) Indian, and (D) Pacific Ocean. Same colorbar is used for panels (A–D) as shown below them. The contour interval is 10.0 mol/m². (E) Annual zonally-averaged O₂C for the world ocean and the three major basins. Monthly time series estimate of the O₂C between 30° and 70° latitude for the (F) World, (G) Atlantic, (H) Indian and (I) Pacific Ocean, respectively. The shadow areas in (F–I) represent the data variability (±1 SD). Note O₂C in this figure is area-averaged to compare the values of the world ocean with the three basins.

Oxygen content shows the total O₂ mass variations in the surface mixed layer, which is highly correlated with the mixed layer thickness (depth). The average O₂ concentration in the mixed layer, O_2m, is a variable that less depends on the mixed layer depth. In Figure 3, we show the zonally-averaged mean O₂ concentration in the ML for the global ocean and the three major basins. Interestingly, the seasonal variations and spatial distributions of O_2m are very similar for the global ocean and the three major basins (Figures 3A–D), even though their respective O₂C distributions (Figure 2) are quite different. The annual averaged O_2m for the global and the three basins are also very similar and mostly a function of latitude (Figure 3E). That is because the water in the mixed layer is saturated or near-saturated (O₂ saturation >90% in the ML) for the global ocean and the three basins (Figure S1 in the Supplemental Material) throughout the year in the mixed layer and the mean O₂ concentration becomes mainly a function of the O₂ solubility of the water, which is lowest at tropical regions and increases towards polar regions (Keeling and Garcia, 2002). Although the O₂ solubility is a function of temperature and salinity, temperature is the main driver for the latitudinal O₂ solubility increase towards polar regions because temperature field in the mixed layer is also a function of latitudes (Boyer et al., 2018; Garcia et al., 2018; Locarnini et al., 2018). The maximum O_2m is typically observed in May in the NH and in December or January in the SH (Figures 3F–I) when the oxygen saturation reaches 100% (Figure S1). The minimum O_2m usually happens in September in the NH (Figures 2F, G, I). The contour lines for the O_2m in the SH are flatter between 30°S and 60°S than those in the NH and the seasonal variations are not as obvious as those for the NH. O_2m minimum is observed between September and November south of 60°S (Figure 3), when and where the oxygen saturation is lowest (Figure S1). The time series in Figures 3F-I also showed similarities regarding the seasonal variations in the global and basin averaged O_2m.

FIGURE 3

Figure 3 Climatology of mean O₂ concentration in the surface ML (mol/m³). Zonally averaged monthly mean oxygen concentration in the surface mixed layer for the (A) World, (B) Atlantic, (C) Indian, and (D) Pacific Ocean. Same colorbar used for panels (A–D) as shown below them. The contour interval is 0.02 mol/m³. (E) Annual zonally averaged oxygen for the world ocean and the three major basins. Monthly time series estimate of the mean oxygen between 30° and 70° for both northern Hemisphere (solid black) and Southern Hemisphere (dashed black) for the (F) World, (G) Atlantic, (H) Indian and (I) Pacific Ocean, respectively. The shadow areas in (F–I) represent the data variability.

3.3 Mixed layer oxygen budget

In this section, we apply the O₂ mass balance model (Section 2.4) to study the controlling factors of O_2m and O₂C annual cycle for the global ocean. Entrainment, air-sea oxygen flux and “biological” O₂ production are the three major controlling factors that are considered in this simple model to examine the relative importance of each factor in each month. Figures 4A–D show the zonally-averaged monthly values for the four terms in Eq. (2). For all the terms examined, positive means adding oxygen to the mixed layer and negative means losing oxygen from the mixed layer, either to the deeper layer or to the atmosphere. In the case of air-sea flux, negative means outgassing and positive ingassing. The seasonality of the four terms is most significant between 30˚ and 70˚ for both NH and SH (Figure 4). At high and mid-latitudes, the oxygen annual cycle can be separated into four stages based on the four terms in Figure 4. During the spring, when the mixed layer depth shoals (Figure 1), the oxygen content begins to decrease with time from the maximum content values (Figures 2A, F, 4A, E). The entrainment term also decreases and weakens to zero due to the shoaling of the mixed layer (Figures 4B, F). The air-sea O₂ ingassing weakens and a reversal to outgassing occurs at middle to late spring (Figures 4C, G). The biologically-produced oxygen increases rapidly and peaks in late spring (Figures 4D, H). In the summer, the oxygen content change is slightly negative and small (Figure 4A). The outgassing of oxygen is balanced with the oxygen produced by biological activity (Figure 4). In the fall, outgassing transitions to ingassing due to solubility changes (colder = more soluble) and oxygen starts to transfer downward via entrainment (Figures 4B, C). Biological processes change from a net source (production) to a net sink (consumption). During the winter, the entrainment continues transferring a large amount of oxygen into the deepwater. The biological activity also uses up oxygen in the ML. The air-sea ingassing is the primary source of oxygen in the ML at high and mid-latitudes (>30˚ latitude). In the tropical regions between 30˚S and 30˚N, the annual net oxygen content change and entrainment are near zero (Figure 4I) suggesting that the air-sea flux is nearly balanced by the biological activity and other forcing processes that are not considered in Eq. (2).

FIGURE 4

Figure 4 Zonally averaged mixed layer (A) oxygen content change, (B) entrainment, (C) air-sea O₂ flux and (D) residuals (“biological” O₂ production). Same colorbar is used for panels (A–D) as shown below them. Monthly time series estimate between 30° and 70° for (E) oxygen content change, (F) entrainment, (G) air-sea O₂ flux and (H) “biological” O₂ production, respectively. (I) Annual zonally-averaged net O₂ flux (change) for the terms shown in (A–D). The shadow areas in (E–H) represent the data variability (±1 SD).

In the high and mid-latitudes of the NH, the period of positive oxygen content changes lasts for about 8 months (from September until April, Figure 4A), which yields a net annual oxygen increase in the ML between 25°N and 75°N (Figure 4I). The annual oxygen content change between 25°S and 25°N is about zero. In the SH, the annual oxygen content change is positive between 25°S and 62°S, but negative between 62°S and 75°S. The annual entrainment flux is always negative and only significant in the high and mid-latitudes (Figure 4I). The annual net air-sea O₂ flux is predominately ingassing in the high latitudes and outgassing in the mid and low latitudes. The annual net biological O₂ production in the ML is positive for most regions of the global ocean except for the high latitudes near Antarctic, suggesting biological activity contributes net O₂ to the ML on a global scale.

To better demonstrate the relative importance of the three controlling factors on the right-hand side of Eq. (2), in Figure 5, we replot the monthly time series shown in Figures 4E–H using histograms. Figures 5A, C display the actual flux values in mole m^-2 yr^-1 in the NH and SH, respectively; while Figures 5B, D show the relative percentage (importance) of the magnitude of the three factors by scaling their total magnitude to 100% in each month. We also superimpose the mixed layer depth on the histograms in Figures 5A, C to show the seasonality of the MLD. In the winter, air-sea flux is the dominant factor, controlling 50%-70% of the oxygen content change in both NH and SH. Biological activity becomes a dominant controlling factor in early spring (50-80%). Between late spring and early fall, the air-sea flux plays a slightly larger role (50-60%) than the biological activity (40-50%). The entrainment effect is only strong from late fall to early spring (contributing to 30%-40% of the O₂C change, Figure 5). Air-sea flux is positive (ingassing) in the winter and negative (outgassing) in the summer in most regions (Figure 4C). In most months, it controls 40%-60% of the relative O₂C change (Figure 5). Biological activity produces oxygen in spring and summer and consumes oxygen in fall and winter in the ML based on the J_o2 estimate (Figure 4D). Biological activity plays a relatively larger role in the spring/summer (40-50%) than that in the fall/winter (10-30%, Figure 5).

FIGURE 5

Figure 5 (A) Monthly bar-chart of the four terms in Eq. (2) for the Northern Hemisphere (NH) between 30°N and 70°N. This is a replot of the monthly time series shown in Figures 4E–H. Superimposed black line is the average mixed layer depth between 30°N and 70°N to show the seasonality of the MLD with the four terms. (B) Relative percentage (RP) of the magnitude (absolute value) of the three terms on the right-hand side of Eq. (2) to show the relative importance of each term in controlling O₂C change in each month in the NH between 30°N and 70°N. Blue: entrainment, red: air-sea flux, and green: “biological” oxygen production. (C) Monthly bar-chart of the four terms in Eq. (2) for the Southern Hemisphere (SH) between 30°S and 70°S. Similar notations are used as in (A). (D) RP of the magnitude of the three terms in Eq. (2) between 30°S and 70°S. Similar notations are used as in (B).

Based on the above analysis, the oxygen content and O_2m annual cycle in the surface mixed layer is controlled by multiple complicated processes that non-linearly correlated with the mixed layer depth. Some of the processes are directly related to the mixed layer depth annual cycle (e.g. entrainment). However, some of the processes are more determined by other factors (e.g. air-sea flux and biological activity), which only indirectly relate to the mixed layer depth changes. The thickness of the mixed layer is determined by wind forcing and heating from the atmosphere. Wind forcing is also a key parameter determining the air-sea oxygen flux at sea surface (see Section 2.4.2). Heating from the atmosphere affects the water temperature, which further affects the oxygen solubility and mixed layer depth. Variations in mixed layer affect the rate of air-sea oxygen flux between the atmosphere and deeper ocean, the capability of the ocean to store oxygen and nutrient supplies to support the growth of phytoplankton. All these processes interact each other and build up complex linkages between mixed layer depth and oxygen concentration/content annual cycle in the upper ocean. The simple mass balance model effectively demonstrates the relative importance of the three essential factors in different seasons for the upper ocean oxygen dynamics in the global ocean. This is also valuable for improving our understanding of overall biology activity in the surface mixed layer. In the following section, we will discuss more about biologically-produced oxygen and relate it to the net community production in the mixed layer for the global ocean.

3.4 Annual net community production for the global ocean

The seasonality of J_O2 follows our traditional understanding of biologically-produced oxygen patterns (producing oxygen in spring/summer and consuming oxygen in fall/winter). However, it is difficult to directly validate the J_O2 estimates due to the lack of biological O₂ production observations in the mixed layer for the global ocean. We convert J_O2 to net community production (NCP) following Craig and Hayward (1987); Emerson (1987); Spitzer and Jenkins (1989); Emerson et al. (2008) and Yang et al. (2017) and compare our results with previous NCP estimates. Assuming a ΔO₂:ΔC ratio of 1.45 (Hedges et al., 2002; Emerson and Stump, 2010; Bushinsky and Emerson, 2015), J_O2 is converted to NCP and shown in Figure 6. Our estimates in Figure 6 provide a global view of the seasonal variations of the NCP for the world ocean based on the oxygen mass balance. Our estimates of zonal annual net community production (ANCP, Figure 6B) range from approximately 0.5 mol C m^-2 year^-1 to ~6.0 mol C m^-2 year^-1 except for south of 60˚S, which can be less than -5.0 mol C m^-2 year^-1. NCP in mid and high latitudes is largest in the spring and smallest in the winter (Figures 6A, C), consistent with the typical seasonal biological activity in the ocean. The global average ANCP is 3.06±1.61 mol C m^-2 year^-1 (mean±SD) between 60˚S and 60˚N. The average ANCP in the NH between 0 and 60˚N is 3.69±1.38 mol C m^-2 year^-1, which is much larger than the average value (2.41 ±1.57 mol C m^-2 year^-1) in the SH (0-60 ˚S). Similar hemispheric asymmetry has been reported in Quay et al. (2009) in the Pacific with subtropical ANCP in the north Pacific (2.6±0.6 mol C m^-2 year^-1) being twice that in the South Pacific (1.3±0.4 mol C m^-2 year^-1). Our results suggest the ANCP hemispheric asymmetry may be a global phenomenon. The causes of asymmetric ANCP distribution between NH and SH are unclear and could inspire further investigation. Our global average ANCP estimates are comparable to previous studies at fixed stations or in regional seas. ANCP at different long-term fixed stations are estimated to be from 1.9-3.8 mol C m^-2 yr^-1 from 10+ year time series (Keeling et al., 2004; Emerson, 2014). Stanley et al. (2010) calculated the mean NCP in the equatorial Pacific to be 1.2~1.5 mol C m^-2 yr^-1 based on data from a 2-month summer cruise. Yang et al. (2017) estimated the ANCP in the subtropical Pacific Ocean to be 0~2.4 mol C m^-2 yr^-1 using Argo oxygen data. ANCP estimates based on models range from 1.0 to 5.7 mol C m^-2 yr^-1, while satellite-based ANCP estimates range from 1.5 to 3.8 mol C m^-2 yr^-1 (Emerson, 2014). Haskell et al. (2020) reported similar total carbon export out of the euphotic zone (ANCP) of 3.2±0.6 to 4.5±0.5 mol C m^-2 yr^-1. Note, unlike using a variable mixed layer depth as in this study, in the previous ANCP estimates, fixed integral depth (e.g. deepest annual MLD, euphotic zone) is often used in the ANCP estimates, which makes it difficult to conduct a direct comparison among studies because different integral depths are used. Moreover, different assumptions and equations were used in different estimate methods, which makes the comparison further complicated. Nonetheless, all these estimates are for the ANCP in the upper ocean and a rough comparison can be still conducted by comparing the range and trend in different estimates. Large discrepancies might be seen when comparing ANCP from different estimate methods. Keep in mind, the ANCP estimated in this study is the ANCP within the mixed layer of the global ocean. A portion of the annual productivity between summer MLD and euphotic depth is not included in our ANCP estimate when the MLD is shallower than the euphotic depth.

FIGURE 6

Figure 6 (A) Zonally averaged mixed layer net community production (NCP) determined by the O₂ mass balance model. (B) Zonally-averaged annual net community production (ANCP). (C) Monthly time series estimate of NCP between 30° and 70° for both NH and SH. The shadow areas in (B, C) represent the data variability (±1 SD).

From our mass balance model, the global ANCP (60˚S-60˚N) in the mixed layer is estimated at 859.7±73.8 Tmol C yr^-1, which is slightly larger than the estimate (781±393 Tmol C yr^-1) from a satellite-based model in the euphotic zone (Westberry et al., 2012). ANCP is difficult to quantify on a global scale and estimates of ANCP (or sometimes carbon export production) based on different methods vary considerably, from 250 to more than 2800 Tmol C yr^-1 (Sherr and Sherr, 1996; Laws et al., 2000; del Giorgio and Duarte, 2002; Najjar et al., 2007; Emerson, 2014). Our global ANCP estimate is within the range, but close to the lower end compared to some of the previous estimates. This is likely due to the fact that a larger criterion (0.125 kg m^-3) is used in our MLD definition. Integrating over a larger mixed layer thickness may include more community respiration in our calculation than some of the previous studies using 0.03 kg m^-3 as their MLD criterion or surface satellite-based observations.

In the SH, between 30˚S and 60˚S, our ANCP estimates agree with previous calculations in Johnson et al. (2017) and Arteaga et al. (2019) in demonstrating a low value in subtropical region (30˚S-40˚S) and an increase in ANCP in the Polar Front Zone and Polar Arctic Zone (40˚S-55˚S, Figure 6B). Our negative J_O2 and ANCP estimate near Antarctic with latitudes greater than 62˚S (Southern Seasonal Ice Zone) is consistent with the estimate in Trucco-Pignata et al. (2022) in the southeast Pacific, but quite different than that in Johnson et al. (2017), Arteaga et al. (2019) and Yang et al. (2021). Although the lowest, but still positive ANCP (1-2 mol C m^-2 yr^-1) was calculated by Johnson et al. (2017) and Arteaga et al. (2019) between 60˚S and 70˚S. Omitting diapycnal eddy diffusion might be one of the reasons why a negative value is calculated in our ANCP in the Seasonal Ice Zone, south of 62˚S. Yang et al. (2021) found the diapycnal eddy diffusion played an important role in determining the net annual mixed layer dissolved inorganic carbon on the West Antarctic Peninsula Shelf based on annual carbon budget from a subsurface mooring. Another reason of the negative ANCP might be associated with the large amplitude ingassing air-sea O₂ flux (Figure 4). The Southern Ocean is an O₂ unsaturated area throughout the year (Figure S1), which makes the diffusive O₂ flux ingassing year-long (Ito et al., 2004). The magnitude of the wind is also relatively large throughout the year in that region, which helps the bubble induced O₂ into the ocean. More investigation and observational data are needed to further examine the ANCP in the Southern Seasonal Ice Zone.

Despite the overall consistency of the seasonal variations in oxygen content and ANCP with former studies, there are a number of improvements that are needed and there are limitations in this simple oxygen mass balance model. Our air-sea oxygen flux is an update to the estimate made by Najjar and Keeling (2000). 20 more years of oxygen data collected after Najjar and Keeling (2000) paper are added in deriving the WOA18 oxygen climatology. The most updated wind fields, sea surface pressure and sea surface skin temperature monthly climatologies are used in our calculations, which should be more accurate benefiting from the development in satellite observations in recent 20 years. Furthermore, new bubble models that recently developed by Liang et al. (2013) and further improved by Emerson et al. (2019) are used to estimate the bubble-mediated oxygen ingassing, which is a significant improvement in estimation of air-sea gas flux. However, although Garcia et al. (2010) showed the standard error of data used to derive the WOA climatology in the majority of the open ocean could reach to 1-2%Sat, the accuracy of WOA18 oxygen climatology may still not high enough to precisely capture the fine-scale air-sea gas flux (Takeshita et al., 2013). The global oxygen observations to date are still too sparse, especially during winter time and in high latitudes. More observations are needed to improve the accuracy of WOA oxygen climatology.

In this simple model, both diapycnal and horizontal diffusion effects are not considered which might cause uncertainties and errors in the J_O2 calculation. Larger turbulent diapycnal mixing is estimated during fall/winter than that in spring/summer, suggesting seasonal variations might exist in turbulent diffusivity based on heat and salt budgets in the North Pacific (Cronin et al., 2015). This also suggests that the errors due to ignorance of diapycnal mixing might be larger in fall/winter. Ignoring the horizontal advection and diffusion term in the model might also cause large errors/uncertainties in the frontal regions (e.g. Kuroshio Extension, Bushinsky and Emerson, 2018), which should be considered when studying the finer-scale studies in addition to the basin-wide and global-scale analyses. The uncertainties and errors in the O₂ climatology, wind and current fields may also affect the J_O2 estimates, especially in the regions with sparse O₂ observations. A full oxygen mass budget may be needed in the future to evaluate the advection, diffusion and other effects (e.g. as in Bushinsky and Emerson, 2015; Bushinsky and Emerson, 2018) for the global ocean.

4 Summary

The spatial and seasonal evolutions of O₂C and O_2m in the surface ML are determined for the global ocean and the three major basins. An O₂ mass balance model is utilized to quantify the entrainment, air-sea flux and biological activity in controlling the seasonal O₂ content/concentration change. Based on the mass balance, from mid-fall until early spring in mid and high latitudes, the ML oxygen change is dominated by ingassing from air-sea flux and oxygen loss due to entrainment and biological consumption. The oxygen content change reaches a maximum in the spring due to increased biological oxygen production, while entrainment vanishes and air-sea exchange transitions from ingassing to outgassing. The oxygen content decreases with time due to the strong outgassing in the summer, which overcomes the oxygen produced by biological activity. The entrainment effect is insignificant in the summer. The model estimated mean ANCP ranges from ~0.5 mol C m^-2 year^-1 to ~6.0 mol C m^-2 year^-1 with an average value of 3.06±1.61 mol C m^-2 year^-1 between 60˚S and 60˚N latitude, comparable to previous estimates. The global ANCP is estimated at approximately 863.7±73.8 Tmol C yr^-1, which is in good agreement with the estimate (781±393 Tmol C yr^-1) based on the photosynthesis/respiration ratio and satellite-based observations (Westberry et al., 2012). The analyses presented here provide a global view at the factors controlling the oxygen concentration in the surface mixed layer and also provide a global estimate of the annual net community production based on the oxygen mass balance.

Data availability statement

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

Author contributions

ZW and HG contributed to conception and design of the study. ZW conducted mass balance model development and data analysis. JR developed the mixed layer depth climatology. TB and HG led the WOD and WOA management and updates. ZW wrote the first draft of the manuscript. JC contributed to funding acquisition. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by NOAA NCEI and a NOAA grant 363541-191001-021000 (Northern Gulf Institute) at Mississippi State University.

Acknowledgments

The authors would like to thank the scientists and data providers for continuing to share their data, data managers and our NCEI colleagues for stewarding the data for long-term preservation, and the Ocean Climate Laboratory colleagues for continuing to enhance and maintain the WOD and WOA. Without this combined effort, data-based studies like this would not be possible. We thank Y. Lau for helping on air-sea oxygen flux calculation. We thank S. Cross, C. Paver, C. Bouchard, K. Larsen, H. Zhang and other NCEI colleagues for valuable conversations and comments.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022.1001095/full#supplementary-material

References

Andreev A. G., Baturina V. I. (2006). Impacts of tides and atmospheric forcing variability on dissolved oxygen in the subarctic north pacific. J. Geophys. Res.: Oceans 111 (C07S10), 1–11. doi: 10.1029/2005JC003103

ORIGINAL RESEARCH article

Controlling factors of the climatological annual cycle of the surface mixed layer oxygen content: A global view

1 Introduction

2 Data and methods

2.1 Climatological oxygen fields and WOA18

2.2 Climatological mixed layer depth

2.3 Oxygen content calculation

2.4 Oxygen mass balance model

2.4.1 Entrainment FE

2.4.2 Air-sea Oxygen Flux FA-W

3 Results and discussion

3.1 Zonal and seasonal variation of MLD

3.2 Zonal and seasonal variation of O2C and O2m in the ML

3.3 Mixed layer oxygen budget

3.4 Annual net community production for the global ocean

4 Summary

Data availability statement

Author contributions

Funding

Acknowledgments

Conflict of interest

Publisher’s note

Supplementary material

References

2.4.1 Entrainment F_E

2.4.2 Air-sea Oxygen Flux F_A-W

3.2 Zonal and seasonal variation of O₂C and O_2m in the ML