Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences
You have accessResearch article

Biogeochemical modelling of dissolved oxygen in a changing ocean

Oliver Andrews

Oliver Andrews

Tyndall Centre for Climate Change Research, School of Environmental Sciences, University of East Anglia, Norwich NR4 7TJ, UK

[email protected]

Google Scholar

Find this author on PubMed

,
Erik Buitenhuis

Erik Buitenhuis

Tyndall Centre for Climate Change Research, School of Environmental Sciences, University of East Anglia, Norwich NR4 7TJ, UK

Google Scholar

Find this author on PubMed

,
Corinne Le Quéré

Corinne Le Quéré

Tyndall Centre for Climate Change Research, School of Environmental Sciences, University of East Anglia, Norwich NR4 7TJ, UK

Google Scholar

Find this author on PubMed

and
Parvadha Suntharalingam

Parvadha Suntharalingam

School of Environmental Sciences, University of East Anglia, Norwich NR4 7TJ, UK

Google Scholar

Find this author on PubMed

    Abstract

    Secular decreases in dissolved oxygen concentration have been observed within the tropical oxygen minimum zones (OMZs) and at mid- to high latitudes over the last approximately 50 years. Earth system model projections indicate that a reduction in the oxygen inventory of the global ocean, termed ocean deoxygenation, is a likely consequence of on-going anthropogenic warming. Current models are, however, unable to consistently reproduce the observed trends and variability of recent decades, particularly within the established tropical OMZs. Here, we conduct a series of targeted hindcast model simulations using a state-of-the-art global ocean biogeochemistry model in order to explore and review biases in model distributions of oceanic oxygen. We show that the largest magnitude of uncertainty is entrained into ocean oxygen response patterns due to model parametrization of pCO2-sensitive C : N ratios in carbon fixation and imposed atmospheric forcing data. Inclusion of a pCO2-sensitive C : N ratio drives historical oxygen depletion within the ocean interior due to increased organic carbon export and subsequent remineralization. Atmospheric forcing is shown to influence simulated interannual variability in ocean oxygen, particularly due to differences in imposed variability of wind stress and heat fluxes.

    This article is part of the themed issue ‘Ocean ventilation and deoxygenation in a warming world’.

    1. Introduction

    The representation of realistic oxygen fields is an important and active area for ocean biogeochemistry model development (e.g. the Fast Met Office UK Universities Simulator [FAMOUS] [1]) owing to the ecological importance of dissolved oxygen concentration (hereafter [O2]) and its sensitivity to climatic perturbation on interannual to millennial time scales [2]. Moreover, despite occupying less than 10% of the ocean by volume, low-oxygen waters could account for up to 50% of open ocean N2O production as a consequence of denitrification and increased N2O yields from nitrification under low-[O2] conditions [3]. Earth system models (ESMs) consistently project a reduction in the oxygen content of the global ocean in response to future anthropogenic climate change [4,5], but with significant uncertainty in the spatial pattern (or ‘fingerprint’) and magnitude [6,7]. Therefore, the ability of current models to reproduce the observed dynamics and distribution of ocean [O2] needs to be scrutinized to understand uncertainties and improve projections [8,9].

    A number of persistent model biases have been identified in the [O2] fields of ocean biogeochemistry models of differing complexity [9,10] integrated for the historical period. Most evident is the inability of current models to reproduce the observed distribution and variability of [O2] at low latitudes, particularly within oxygen minimum zones (OMZs) such as the eastern tropical Pacific [11]. For example, Ocean Biogeochemistry General Circulation Models (OBGCMs) generally simulate historical [O2] increases within the tropical thermocline (e.g. [12,13]) and an overall contraction of suboxic waters [7] in response to ocean warming. These model responses are contrary to time-series data compiled for the eastern tropical OMZs [14], which show marked deoxygenation trends and expansion of low-[O2] waters over the last approximately 50 years. Gnanadesikan et al. [15] and Cabré et al. [11] also highlight the limited capacity of current ESMs to reproduce the observed climatological distribution of [O2] in the subsurface. This systematic model bias is entrained into the biogeochemically significant [16] volume censuses of low-[O2] waters, with models variously both underestimating (e.g. HadGEM2-ES [17]) and overestimating (e.g. MPI-ESM-LR [18]; GFDL-ESM2.1 [19]) hypoxic and suboxic extent within the ocean interior. A range of dynamical and biogeochemical model deficiencies have been invoked to explain the divergence between observed and simulated [O2] at low latitudes (e.g. [1,9]). Chief among them is the inability of coarse resolution models to explicitly simulate the mesoscale structures which resupply O2 into the eastern tropical OMZs (e.g. [20]). Eddy-resolving spatial resolution has been shown to improve representation of [O2] in the Arabian Sea [21] and eastern tropical Atlantic [10] OMZs via more realistic transport processes in the physical model. However, below eddy-resolving scales, improvements in the representation of [O2] fields with increasing spatial resolution have been shown to be small (e.g. MPI-ESM [18]). Uncertainties surrounding the parametrization of lateral [15,22] and vertical [23] diffusion have also been shown to place important constraints on the extent of low-[O2] waters in coarse-resolution ESMs. In addition, unforced variability in the equatorial trade winds has been shown to influence the simulated evolution of tropical OMZs on multidecadal time scales [24], further complicating model–data agreement in coupled models.

    The representation of biological processes and their response to physical and geochemical forcing has been suggested as an important limitation on O2 dynamics in ocean biogeochemistry models (e.g. [25,26]). For example, Matear & Lenton [27] show, using an ensemble of ESM experiments, that the interactive effects of climate change and ocean acidification could drive significant alteration to biogeochemical cycles, suggesting that ocean warming should not be considered in isolation when considering future ocean deoxygenation. The potential impacts of ocean acidification on primary and export production are particularly important for oxygen cycling. Mesocosm experiments carried out using natural plankton communities [28,29] suggest that the C : N uptake ratio in photosynthetic carbon fixation increases under elevated pCO2 due to dissolved inorganic carbon (DIC) overconsumption, causing the composition of exported organic material to deviate from classical Redfield stoichiometry (e.g. C : N = 6.6 [30]). Specifically, Riebesell et al. [28] report higher C : N drawdown ratios in diatom and coccolithophore-dominated mesocosm enclosures exposed to increased partial pressures of CO2 (pCO2 = 700 µatm, C : N = 7.1; pCO2 = 1050 µatm, C : N = 8.0), while N : P ratios remain unchanged from Redfield proportions. The effects of stoichiometric plasticity in marine ecosystems remain largely unaccounted for in current ocean biogeochemistry models [31], which generally rely on fixed elemental ratios in the formation of organic material. However, studies using configurations of the UVic [32] and PISCES [33] models, which implement a pCO2-sensitive C : N ratio in primary production, suggest that stoichiometric effects in response to ocean acidification could have a major impact on biogeochemical cycles. Model experiments with variable C : N stoichiometry simulate increases in cumulative carbon export of between 70 and 100 Pg C by 2100 in response to future CO2 forcing. Moreover, the export of (relatively) more carbonaceous organic material from the surface in response to DIC overconsumption driven by acidification also causes biological oxygen demand to increase. As a result, these model experiments project elevated deoxygenation of the tropical thermocline (greater than 20 µmol l−1) and 36–50% increases in the volume of suboxic waters by 2100 [32,33], compared with no change or a small [O2] increase as simulated by fixed stoichiometry ESMs [6]. Stramma et al. [9] investigate the impact of variable C : N ratios on historical [O2] trends using an intermediate complexity ESM, and find that this effect is small for zonal mean [O2] trends at 300 dbar. However, the impact of accounting for variable stoichiometry on historical [O2] changes in more complex plankton functional type (PFT) ocean biogeochemistry models (e.g. [34,35]) remains uncertain.

    Another potentially important impact of biological processes on [O2] changes relates to the well-established [36] effect of ocean acidification on the saturation state of seawater with respect to calcium carbonate (both calcite [Ωcal] and its less stable polymorph aragonite [Ωarag]). A decrease in Ω under elevated pCO2 has been shown by a number of laboratory and field studies to reduce biogenic calcification rates in some calcareous holoplankton (foraminifera, pteropods and some coccolithophores) and warm-water corals [37]. Extrapolating these results, Heinze [38] projects an approximately 50% decline in global CaCO3 export production by 2250 using an OBGCM which accounts for reduced biocalcification rates in response to ocean acidification (SRES A1B emissions scenario). Beyond direct effects, it has been suggested by Armstrong et al. [39] that changes in the ‘carbonate pump’ may also impact upon the export of organic material, because fluxes of particulate organic carbon (POC) and particulate inorganic carbon (PIC) are highly correlated beneath approximately 1000 m depth (the so-called PIC : POC rain ratio [40,41]). Specifically, export fluxes of dense calcareous (CaCO3) and siliceous (SiO2) biominerals provide ballast which increases the sinking speed and transfer efficiency of POC into the ocean interior. Therefore, reduced export production of CaCO3 in response to ocean acidification could impact upon the efficiency of the organic (soft tissue) biological pump, such that POC remineralizes at shallower depths as mineral ballast fluxes weaken. This effect has been reproduced by OBGCMs that include simple ballasting sub-models [38], with implications for other biogeochemical cycles including oxygen. For example, Hofmann & Schellnhuber [42] demonstrated, using idealized experiments with a [Inline Formula] dependency in calcification rates (RCAL), that shallower POC remineralization depths in response to weakened ballasting exacerbates O2 depletion within established subsurface [O2] minima. Specifically, reduced ballasting alone is shown in prognostic experiments under a SRES A1F1 emissions scenario until 2100 (with emissions declining to 0 at 2200) to reduce [O2] by 20–40 µmol l−1 between approximately 200 and 800 m depth by the year 3000, with the largest decreases within tropical OMZs (e.g. greater than 50 µmol l−1 at 500 m depth in the Arabian Sea). Thus, the effects of ocean acidification on carbon export have been shown to have a major impact on [O2] distributions and forced responses in model projections. However, no studies have so far addressed if these missing processes could reconcile the mismatch between observed and modelled [O2].

    Current ocean biogeochemistry models also generally underestimate temporal variability in O2 on interannual to decadal time scales. For example, Deutsch et al. [43,44] report decadal scale variability in apparent oxygen utilization (AOU) as simulated by an ensemble of ocean-only hindcast model experiments for the North Pacific region to be underestimated by a factor of approximately three relative to repeat hydrographic section data. Similarly, interannual variability in global atmospheric potential oxygen fluxes have been shown to be underestimated by a factor of approximately two to four in hindcast ocean biogeochemistry models relative to those estimated using global atmospheric transport inversions [26]. These results are also consistent with the model-data comparison of [8], where optimal detection methods applied to Coupled Model Intercomparison Project Phase 5 (CMIP5) ESMs show that model [O2] responses need to be scaled up by a factor of approximately two to four to match observed changes. Recent work using hindcast models [4547] demonstrates the sensitivity of simulated trends and variability in the ocean carbon cycle to imposed atmospheric forcing. For instance, Ishi et al. [47] find significant differences in the interannual variability of CO2 outgassing fluxes from the tropical Pacific using ocean-only OBGCM experiments forced with the JPL CCMP Ocean Surface Wind Product [48] compared with NCEP/NCAR reanalysis data [49]. Additionally, model representation of large-scale physical transport processes, such as the Atlantic Meridional Overturning Circulation (AMOC), has been shown to be sensitive to the choice of imposed atmospheric forcing. For example, Stepanov & Haines [50] show model-simulated long-period variability in AMOC transport at 26.5°N to be significantly different between hindcast experiments which employ ECMWF ERAInterim reanalysis data [51] compared to those which calculate bulk fluxes using the Drakkar Forcing Set 3 (DFS3; [52]) blended meteorological and satellite forcing product. As a result, while it must be acknowledged that natural variability is also generated internally to the ocean system, ocean-only model configurations allow considerable ‘exogenous’ variability to be directly related to the imposed atmospheric forcing. Therefore, the use of high-frequency, high-quality atmospheric data to calculate turbulent fluxes of heat, freshwater and momentum in hindcast models (cf. [52,53]) could provide a mechanism for improving interannual to decadal variability in simulated O2.

    The present study explores the source of systematic biases in model representation of [O2], with a range of physical and biogeochemical perturbation experiments conducted using a state-of-the-art global ocean biogeochemistry model. It focuses on quantifying the impact of (i) ocean acidification (variable C : N stoichiometry and mineral ballasting) and (ii) imposed atmospheric forcing on the spatio-temporal distribution of [O2] over the last approximately 50 years.

    2. Model description

    (a) Ocean biogeochemistry model

    PlankTOM10 is a global ocean biogeochemistry model, which describes lower trophic-level ecosystem dynamics explicitly based on PFTs [35,54,55]. It includes 10 PFTs: picophytoplankton, N2-fixers, coccolithophores, mixed phytoplankton, diatoms, colonial Phaeocystis, bacteria, protozooplankton, mesozooplankton and macrozooplankton. The present model version, fully documented in [35], comprises 39 biogeochemical tracers and simulates the full marine cycles of carbon (C), phosphorous (P), oxygen (O2), silicon (Si), along with simplified cycles of nitrogen (N) and iron (Fe). Growth of PFTs is co-limited by temperature, light, macronutrients (N, P and Si) and Fe. PlankTOM10 includes three detrital pools (semi-labile dissolved organic material, and small and large POC), with a fixed Redfield stoichiometry (172O : 122C : 16N : 1P [56]) in the formation and remineralization of organic material for the standard model configuration. Ratios for Fe : C, Chl : C and Si : C (for diatoms) of organic material are variable, calculated by the model based on PFT and abiotic factors [34]. N pools are also subject to denitrification and N2-fixation processes. PlankTOM10 also includes mineral ballasting, whereby the sinking speed of large POC increases as a function of opal (SiO2) and calcite (CaCO3) content. This parametrization applies the direct measurements of mineral ballasting by opal and calcite in copepod faecal pellets [57] to the drag equations [58]. The drag equations are solved offline by iteration in order to calculate gravitational sinking speeds of large POC (vsink) over the range of particle densities in the model. These pairs of particle densities and sinking speeds could be well represented by the following equation:

    Display Formula
    2.1
    where a = 0.0303, b = 0.6923, ρsw is the density of seawater and ρpar is the density of the particle, calculated by
    Display Formula
    2.2
    where LPOC is large particulate organic carbon, CAL is sinking calcite, SIL is sinking opal, 240 is wet weight mol−1 POC, 100 is the molar mass of calcite, 60 is the molar mass of opal and ρLPOC = 1.08 kg l−1, ρCAL= 1.34 kg l−1, ρSIL = 1.2 kg l−1 (calculated based on the data of [57]). Small POC is set to sink at a constant rate of 3 m d−1, whereas LPOC has a maximum numerically stable sinking speed set to 150 m d−1.

    Biogeochemical fields are initialized from observations of the Global Ocean Data Analysis Project (GLODAP; [59]) data for DIC and alkalinity, and World Ocean Atlas 2005 for O2 [60] and nutrients [61]. Biological variables are initialized as in [62], with initial concentrations available at http://opendap.uea.ac.uk:8080/opendap/greenocean/Restart/. As described in [35], PFTs were not impacted by initialization and equilibrated within 3 years of model integration. The model is forced with atmospheric CO2 concentration at each timestep [63]. Dust fluxes are interpolated to daily values from the monthly fields of [64], providing surface Fe and Si inputs.

    (b) Physical model

    PlankTOM10 is embedded within the Nucleus for European Modelling of the Ocean v. 3.1 (NEMOv3.1) physical model, which comprises the primitive equation Océan Parallélisé v.9 (OPA9; [65]) ocean GCM coupled to the Louvain–la–Neuve Ice v. 2 (LIM2; [66]) dynamic-thermodynamic sea ice model. This study employs a global configuration (ORCA2; [67]) of NEMOv.3.1 where the model is discretized on a tripolar curvilinear grid with a zonal resolution of 2° and a meridional resolution of 2° × cos(latitude) increasing to approximately 0.5° at the equator and towards the poles. In the vertical, ORCA2 has 30 z-levels with a maximum resolution of 10 m for the upper 100 m, decreasing to approximately 500 m at 5 km depth. Vertical mixing is calculated using a turbulent kinetic energy model [68], with subgrid-scale eddy-induced mixing processes accounted for using the parametrization of [69]. Active tracers are initialized from World Ocean Atlas 2005 temperature [70] and salinity [71] observations. A range of atmospheric forcing data has been used to derive surface fluxes of momentum, heat and freshwater as boundary conditions to the hindcast model, as described in §3b.

    3. Simulation set-up

    As outlined in table 1, ocean biogeochemistry (§3a) and atmospheric forcing (§3b) sensitivity experiments are carried out using the PlankTOM10-NEMOv.3.1 model to explore and review impacts on simulated historical distributions of [O2]. To isolate the impact of perturbations, model fields are generally presented relative to a baseline reference configuration. This approach is consistent with other model studies which evaluate the impact of incorporating new processes [72] or anthropogenic impacts [73] via comparison with an unperturbed experiment. As a result, non-equilibrium artefacts and other dynamical processes are removed, such that this analysis focuses on understanding the impact of these imposed changes on hindcast biogeochemical tracers. This work therefore presents a synthesis of sensitivity experiments and differs from other studies (e.g. [9]) focusing on model-data comparison of absolute changes in [O2] over the historical period. This established approach (e.g. [46]) allows for representation of trends and variability in a ‘process-level’ model study to be improved due to the relatively short hindcast period over which observed vertical gradients in [O2] are preserved.

    Table 1.Summary of the model configurations used in each ocean-only model experiment.

    model experiment time period atmospheric forcing pCO2-sensitive C : N ratio pCO2-sensitive RCAL
    REF 1948–2013 NCEP/NCAR no no
    STO10 1948–2013 NCEP/NCAR yes no
    BAL10 1948–2013 NCEP/NCAR no yes
    CORE2 1948–2007 COREv2-IAF no no
    DFS4 1958–2006 DFS4.3 no no
    IPSL 1948–2005 IPSL-CM5A-LR no no

    Results from ocean acidification experiments are presented averaged over the period 2003–2013 (§4a), allowing for the largest adjustments to initial conditions to occur over the first approximately 50 model years. De-trended results from the ensemble of atmospheric forcing experiments are presented over their common time period (1958–2005) to avoid sampling biases associated with differing temporal ranges (§4b).

    (a) Ocean acidification experiments

    Two biogeochemistry perturbation experiments are conducted over the historical period to investigate the impact of ocean acidification on the spatio-temporal distribution of [O2], alongside a baseline (REF) experiment following the standard model configuration. These comprise explicit representation within PlankTOM10 of (i) a pCO2-sensitive C : N ratio in primary production (STO10) and (ii) a pCO2-sensitive calcification rate (RCAL) and associated impacts on the PIC : POC rain ratio (BAL10). Ocean biogeochemistry experiments are conducted using a common atmospheric forcing (REF; §3b).

    For the STO10 experiment, a change in the C : N ratio of organic matter via increased photosynthetic carbon fixation in response to ocean acidification is parametrized using the results of mesocosm experiments carried out on natural plankton communities under elevated pCO2 [28]. A non-dimensional ‘CO2-sensitivity’ factor is used to provide a relationship between pCO2 model forcing and the C : N ratio in simulated organic carbon production and remineralization. Following the prognostic model study of [32], a linear relationship between imposed pCO2 and C : N is derived based on the results of [28] for measured C : N values provided at pCO2 values of 350, 700 and 1050 µatm. Also following the study of Oschlies et al. [32], this CO2-sensitivity factor is re-scaled back to an assumed pre-industrial pCO2 of 280 ppm which increases the C : N ratio for a given pCO2; therefore, our estimate provides an upper bound for this overconsumption effect over the historical period. In STO10, this factor is multiplied by the rate of net organic carbon production and used to calculate a variable (pCO2-sensitive) O : C ratio in organic material. Variable C : N ratios in organic material are linked to imposed atmospheric pCO2 rather than surface water pCO2, enabling direct comparison with a fixed C : N ratio REF experiment. In support of this, Oschlies et al. [32] show that simulated changes in [O2] over longer (centennial) time scales are insensitive to this distinction.

    For the BAL10 experiment, the impact of a pCO2-sensitive biogenic calcification rate (RCAL) on marine biogeochemical cycles is parametrized using results from laboratory manipulations carried out with coccolithophore species, in which the PIC : POC rain ratio is measured to decrease under elevated [CO2(aq)] [74]. PlankTOM10 includes explicit representation of coccolithophores as a calcifying PFT with growth rates based on observations. Thus, the experimentally derived parametrization can be applied directly to coccolithophore responses rather than more generically. Following the study of Zondervan et al. [74] and Heinze [38], the CaCO3 : CORG production ratio (Inline Formula) is parametrized in BAL10 as a function of pCO2:

    Display Formula
    3.1
    where A = 4.4 × 10−4 (based on [74]), Inline Formula and Inline Formula is prescribed from observations [63].

    (b) Atmospheric forcing experiments

    Gridded atmospheric data are used to prescribe surface boundary conditions to ocean models. In this study, a series of model experiments are conducted using different atmospheric data to investigate the sensitivity of variability in simulated ocean [O2] to imposed forcing. Simulations are carried out using four different atmospheric forcing data products (table 2). These experiments also differ in terms of the bulk formulae used to provide turbulent fluxes of momentum, heat and freshwater to the physical ocean model.

    Table 2.Summary of atmospheric forcing datasets used in model experiments. The temporal frequency of meteorological surface variables is provided in parentheses (di = 6 hourly (diurnal); d = daily, m = monthly).

    exp dataset time period bulk formulation U10, V10, θ, q radsw, radlw precip, snow
    REFa NCEP/NCAR 1948—2013 CLIO NCEP/NCAR (d) n.a.b NCEP/NCAR (d)f
    CORE2 COREv2-IAF 1948–2007 CORE NCEP/NCARc (di) ISCCP-FDc,e (d) GCGCSg (m)
    DFS4 DFS4.3 1958–2006 CORE ERA-40d (di) ISCCP-FDd,e (d) GCGCSg (m)
    IPSL IPSL-CM5A-LR 1948–2005 CORE CMIP5 (d) CMIP5 (d) CMIP5 (m)

    aThe baseline REF atmospheric forcing configuration is used for all ocean biogeochemistry experiments.

    bRadiative fluxes calculated from total cloud cover (tcdc) following [49].

    cBias corrections applied as described in [53,75].

    dBias corrections applied as described in [52].

    eForcing provided as a climatological mean annual cycle for prior to 1984.

    fTotal precipitation rate used without solid fraction (snow).

    gForcing provided as a climatological mean annual cycle for prior to 1979.

    As summarized in table 2, the baseline REF experiment uses a CLIO bulk formulation [76] to calculate surface boundary conditions from daily frequency NCEP/NCAR 10 m air temperature (θ10), 10 m u and v wind components (U10, V10), total precipitation rate (pptn), 10 m specific humidity (q10) and total cloud cover (tcdc). Three other experiments (CORE2, DFS4 and IPSL) are conducted, all of which employ a more recent bulk formulation (CORE; [53,75]), which requires a slightly different set of surface variables, including downwelling shortwave (radsw) and longwave (radlw) radiation and precipitation as a total precipitation rate (precip) and solid fraction (snow).

    CORE2 is forced with Version 2 of the Common Ocean-Ice Reference Experiments Inter-annually Varying Forcing (COREv2-IAF) dataset [53,75]. COREv2-IAF is a hybrid forcing product, which applies corrections to known biases in NCEP/NCAR reanalysis state variables (U10, V10, θ10, q10) and uses satellite-derived radiative flux (ISCCP-FD [77]) and precipitation (merged GCGCS product [75]) estimates so as to limit the imbalance in model heat and freshwater budgets.

    DFS4 is forced with the DRAKKAR Atmospheric Forcing Set v. 4.3 (DFS4.3 [52]) data. Following the COREv2-IAF approach, DFS4.3 is also a hybrid forcing product using satellite radiation (ISCCP-FD) and precipitation (GCGCS). Surface atmospheric state variables are, however, provided from ERA-40 ECMWF reanalysis data [78], with adjustments as described in [52]. ERA-40 is considered a ‘second-generation’ reanalysis product, with improvements in terms of resolution, data assimilation methods and atmospheric models.

    IPSL is forced with model data generated from the output of a CMIP5 ESM. Historical integrations of the IPSL-CM5A-LR ESM [79] were conducted under CMIP5 (1850–2005) with climatic forcings prescribed from observations [80]. Here, output from one ensemble member of the ‘historical’ experiments (r1i1p1) was processed to extract atmospheric data according to the approach of DFS4 and CORE2; however, temporal frequency was limited by the availability of model output. Thus, daily frequency U10, V10, θ2, q2, radsw and radlw were used along with monthly mean precip and snow, with modifications made to the physical model set-up to account for changes in temporal frequency. IPSL-CM5A-LR ESM was selected in order to limit the physical inconsistency with PlankTOM10-NEMOv.3.1, because IPSL-CM5A-LR also uses the NEMO physical model, while also reducing errors associated with spatial interpolation. Other studies have conducted prognostic ocean-only OBGCM experiments using CMIP5 atmospheric data (e.g. [81]); however, the impact of an ESM-derived forcing in hindcast model configurations remains underanalysed.

    4. Results

    (a) Ocean acidification experiments

    For STO10, the effects of ocean acidification on the C : N ratio of organic carbon production drive changes in simulated POC export from the euphotic zone (figure 1a) and DIC distribution (figure 1b,c). Globally, POC export at 100 m is higher by more than 0.2 mol C m−2 yr−1 relative to the fixed C : N ratio REF experiment, corresponding to an approximately 20% area mean increase. The higher POC export is most pronounced (greater than 0.3 mol C m−2 yr−1) within established high-production regions of the global ocean, such as the eastern boundary upwelling system of the equatorial Pacific (figure 1a). Comparable increases in export are also found at mid- to high latitudes (40–60°) within the subpolar North Pacific, North Atlantic and Southern Ocean, where existing high rates of annual primary production associated with the spring blooms are accentuated. By contrast, POC export fluxes show a smaller increase within the mid-latitude oligotrophic gyres; however, this corresponds to a similar relative increase (15–20%) when compared with more productive regions. Coeval historical decreases in DIC at 100 m depth across much of the tropical and mid-latitude ocean (up to 3 µmol l−1) are consistent with carbon overconsumption within the euphotic zone in response to the CO2 fertilization (figure 1b). Lower DIC concentrations ([DIC]) simulated within the euphotic zone are opposed by zonal mean [DIC] increases (greater than 5 µmol l−1) at depth, associated with a strengthened ‘soft tissue’ pump as more exported POC is remineralized (figure 1c). Elevated remineralization-driven alteration to the DIC profile at intermediate depths is associated with high POC sinking speeds and model underestimation of upper ocean bacterial biomass (0–200 m; [35]).

    Figure 1.

    Figure 1. Differences in modelled ocean carbon cycle variables (STO10 simulation minus REF) averaged over 2003–2013. (a) POC export (mol C m−2 yr−1) at 100 m depth, (b) concentration of DIC (μmol l−1) at 100 m depth and (c) zonal mean DIC concentration.

    For BAL10, the global mean ratio of CaCO3 to POC export production (Inline Formula) at 100 m decreases by 4.9% relative to REF, consistent with a reduced rate of biogenic calcification in response to ocean acidification. The largest decreases in Inline Formula occur within the tropical Indian Ocean and eastern tropical Atlantic Ocean, and in mid- to high-latitude regions of the North Atlantic and North Pacific (figure 2a). Reductions in Inline Formula of up to approximately 0.01 are similar in magnitude to those reported by Hofmann & Schellnhuber [42] in [Inline Formula]-sensitive RCAL experiments for 2003–2013. Within the North Pacific, marked decreases in Inline Formula are centred on the North Pacific Current (NPC) region. The eastward flowing NPC is a major transverse surface current which bisects the subtropical and subarctic North Pacific and plays an important role in the resupply of nutrients and oxygen into the interior of the Alaskan gyre [82], where secular [O2] decreases have been observed [83]. Thus, Inline Formula reductions in the NPC could have important downstream implications for biogeochemical cycles in the eastern subpolar North Pacific and California Current region [84].

    Figure 2.

    Figure 2. Differences in modelled ocean carbon cycle variables (BAL10 simulation minus REF) averaged over 2003–2013. (aInline Formula at 100 m depth and (b) gravitational sinking speeds for large POC (vsink, m d−1) between 0 and 2000 m.

    In BAL10, reduced export of CaCO3 mineral ballast from the surface ocean causes a global mean reduction in model-simulated gravitational sinking speeds for large POC (vsink) between 0 and 2000 m depth of 0.2 m d−1 (0.4%). The spatial pattern of vsink reductions is most pronounced (greater than 0.5 m d−1), where Inline Formula decreases are largest (figure 2b). As such, perturbation to the PIC : POC ‘rain ratio’ can be invoked to explain coeval changes in sinking speeds of large POC.

    Changes in the carbon cycle impact upon the spatio-temporal distribution of [O2] in STO10 via changes in the rate of oxygen production and consumption. Mirroring the pattern of DIC concentration changes (figure 1c), zonal mean [O2] in STO10 increases by up to 6 µmol l−1 within the more productive (sub)tropical euphotic zone (figure 3). These near-surface [O2] increases are opposed by marked deoxygenation throughout much of the ocean interior, particularly at mid- to high latitudes, where biological oxygen demand rises as more carbonaceous (higher C : N ratio) organic material is remineralized at depth. Subsurface [O2] depletion reaches a zonal mean maximum of greater than 10 µmol l−1 within intermediate waters of the subpolar North Atlantic (approx. 60°N). The signature of depth-averaged zonal mean [O2] changes scales with latitude, such that the largest [O2] decreases (greater than or equal to 2 µmol l−1) occur poleward of 60° in regions of deep water renewal. As such, the inclusion of acidification effects in organic carbon production could act to augment the historical fingerprint of climate-driven ocean deoxygenation produced by fixed stoichiometry models [6].

    Figure 3.

    Figure 3. Difference in modelled zonal mean [O2] (μmol l−1) compared with simulation REF for STO10 averaged over 2003–2013 (blue colours indicate historical deoxygenation relative to REF).

    The inclusion of a variable C : N ratio in carbon fixation in STO10 also impacts upon the characteristics of simulated low-[O2] waters. Decreases in minimum [O2] values are most pronounced (greater than 8 µmol l−1) within the subpolar North Atlantic and for the Indian Ocean and eastern equatorial Pacific OMZs (not shown). Associated with this intensification of low-[O2] conditions, STO10 simulates an increase in the number of suboxic (+2%, [O2] ≤ 5 µmol l−1) and hypoxic (13%, [O2] ≤ 60 µmol l−1) grid cells, relative to REF. Expansion of low-[O2] waters also impacts upon nitrogen cycling, with the promotion of denitrification processes under suboxic conditions where nitrate is used as an oxidant in the remineralization of organic material. Global area mean denitrification rates increase by 0.27 µmol N m−3 yr−1 (34%, 0–2000 m), associated with oxygen depletion in response to elevated POC export fluxes.

    Compared with STO10, changes in the carbon cycle associated with BAL10 cause almost no change in zonal mean [O2] relative to REF. However, in agreement with the prognostic model results of [42], BAL10 reproduces a small [O2] decrease relative to REF within the ventilated thermocline (100–1000 m, not shown), which is likely to be exacerbated in response to future pCO2 forcing. Aside from imposed pCO2, the more muted magnitude of [O2] change simulated by BAL10 can also be attributed to the relatively small overall global export production which dampens the impact of any acidification-driven reductions in Inline Formula on Inline Formula. However, generally the results of BAL10 suggest that the inclusion of a pCO2-sensitive biocalcification rate in an ocean biogeochemistry model does not impact significantly upon simulated O2 dynamics for the historical period, despite alteration to Inline Formula and sinking speeds of large POC.

    (b) Atmospheric forcing experiments

    Simulated standard deviation in annual mean thermocline (300 m) [O2] is presented as a metric of interannual variability following other model studies which assess the bulk variability properties of upper ocean O2 using either an approximately 300 m depth interval [85] or surface fluxes [25,86]. Variability differs between model experiments (figure 4), with DFS4 exhibiting the largest area mean variability (Inline Formula) when compared with CORE (Inline Formula) and IPSL (Inline Formula). Differences in temporal frequency of input forcing data used for the IPSL CORE experiment (table 2) are found to have only a small impact on simulated Inline Formula (2% change in area mean Inline Formula) in sensitivity analyses carried out with DFS4 at the temporal frequency of IPSL (electronic supplementary material, figure S1). Elevated Inline Formula in the DFS4.3 experiment suggests that ERA-40-derived forcing products generate more exogenous variability in passive tracer fields of ocean-only models when compared with NCEP/NCAR (REF, CORE2)-based atmospheric data. However, a number of other differences between forcing products (table 2) could also contribute to the larger interannual variability in DFS4, for instance, the alteration in surface fluxes caused by referencing of DFS4.3 surface air temperature and specific humidity at 2 m rather than 10 m.

    Figure 4.

    Figure 4. Interannual variability (σ) in modelled annual mean [O2] (μmol l−1) at 300 m depth for a range of model experiments. (a) REF (black in (f)), (b) CORE2 (red), (c) DFS4 (green), (d) IPSL (blue) and (e) ESM (turquoise). [O2] contours are overlain in black for Inline Formula. Zonal mean interannual variability in [O2] at 300 m for all model experiments is presented in (f). A boxcar (low-pass) filter is applied in order to diagnose secular trends in subsurface [O2], with the 10-year running mean being removed at each grid point in order to retain only an estimate of interannual (unforced) variability. The interannual variability is computed over the 1958–2005 time period, common to all forcing products.

    [O2] fields are also taken from biogeochemical output of the IPSL-CM5A-LR CMIP5 ‘historical’ experiment (ESM; figure 4e) used to derive the atmospheric forcing for IPSL. As a result, direct comparison can be drawn between variability in the coupled IPSL-CM5A-LR ESM configuration and the ocean-only PlankTOM10-NEMOv.3.1 experiment conducted using atmospheric fields from this integration. The coupled IPSL-CM5A-LR model (ESM) simulates more variability in thermocline [O2] (area mean Inline Formula) relative to the ensemble mean of all ocean-only atmospheric forcing experiments (Inline Formula). As IPSL and ESM experiments include identical atmospheric forcings along with similar physical and biogeochemical model components (IPSL = PlankTOM10-NEMOv.3.1; ESM  = PISCES-NEMOv3.2 [79]), the residual interannual variability between experiments (Inline Formula) can plausibly be attributed to that which is generated internally to the ocean–atmosphere system under a coupled formulation.

    Simulated interannual variability scales with latitude for all model experiments, such that Inline Formula is most pronounced (σ ≥ 4 µmol l−1) poleward of 40° in both hemispheres. All model experiments exhibit an elevated Inline Formula signal of up to 10 µmol l−1 within the northwestern subpolar gyre of the North Atlantic and subpolar central and western North Pacific (figure 4). These regional patterns have been identified in other O2 modelling studies (e.g. [25,85]) and are associated with the North Atlantic Oscillation (NAO) and Pacific Decadal Oscillation (PDO), respectively, which provide the dominant source of Northern Hemisphere climate variability on interannual to decadal time scales. All experiments converge on a maximum Inline Formula of approximately 6 µmol l−1 at approximately 60°N consistent with variability associated with the dominant climate modes, whereas the signal of interannual variability in the Southern Ocean is less certain—with zonal mean Inline Formula ranging from 1 to 7 µmol l−1 for different atmospheric forcings south of 60°S (figure 4f). Particularly, the CORE2 and IPSL experiments simulate Inline Formula across much of the Southern Ocean, whereas zonal mean Inline Formula exceeds 4 µmol l−1 for REF, DFS4 and ESM (figure 4f).

    Interannual variability in windspeed (σwspd; figure 5a) and near-surface air temperature (σθ; figure 5b) also differ between imposed atmospheric forcing products. DFS4.3 exhibits the largest zonal mean variability in windspeed and near-surface air temperature, consistent with the largest simulated Inline Formula of all ocean-only experiments (figure 4c). However, DFS4.3 windspeed and near-surface air temperature data do not reproduce the meridional structure of σwspd and σθ exhibited by all other forcing fields. Excluding DFS4.3, the largest interforcing divergence in σwspd occurs in the tropics (20°S–20°N) where elevated interannual variability in COREv2-IAF tropical windspeeds (approx. 0.5 m s−1) relative to IPSL-CM5A-LR-derived winds (approx. 0.25 m s−1) can be invoked to explain the more muted tropical variability in thermocline [O2] (Inline Formula) of IPSL compared with CORE2 (figure 4f). By contrast, larger interannual variability in σwspd does not generate a first-order response in Inline Formula at mid- to high latitudes, with, for example, elevated variance in the westerlies over the Southern Ocean for IPSL not being associated with a coeval increase in Inline Formula (figure 4d).

    Figure 5.

    Figure 5. Interannual variability (σ) in de-trended zonal mean (a) windspeed (m s−1) and (b) near-surface air temperature (θ; K) for NCEP/NCAR reanalysis (black), COREv2-IAF (red) and DFS4.3 (green) products, and IPSL-CM5A-LR-derived atmospheric fields (blue). Following [45], windspeed is plotted as a measure of momentum flux rather than wind stress because the latter has a strong dependency on the choice of drag coefficient.

    The largest deviation in zonal mean σθ between atmospheric datasets occurs poleward of 60°S, where COREv2-IAF exhibits lower interannual variability (σθ = ∼0.4 K) when compared with IPSL-CM5A-LR (σθ = ∼0.8 K) and NCEP/NCAR reanalysis (σθ = ∼1.2 K) derived near-surface air temperatures. These σθ differences track the interexperiment divergence in Inline Formula for the Southern Ocean, such that model biases in the simulation of thermocline [O2] variability agree with differences in the magnitude of imposed σθ. Specifically, reduced interannual variability in near-surface air temperature in CORE2 and IPSL can be related to lower Inline Formula relative to REF, which exhibits elevated σθ and, therefore, Inline Formula poleward of 60°S. This modulation of tropical O2 variability by windspeed and extratropical (subpolar) O2 variability by surface heat flux has been reproduced by other forced ocean biogeochemistry models investigating variability in North Atlantic O2 fluxes [86].

    5. Discussion and summary

    To date, much research has focused on understanding the impact of biogeochemical [27,32,33] and physical [6,7] processes in OBGCMs on simulated future climate-driven perturbation to the global oxygen inventory. This work, rather, aims to better constrain the implications of the interrelated physical–biogeochemical drivers for [O2] dynamics over the historical period towards reconciling the well-documented [9] discrepancies between model-simulated and observed [O2] changes. As highlighted by Keeling et al. [2], a ‘critical first step’ in our understanding of oxygen dynamics in a warming world is the development of models at regional and global scales with the capacity to reproduce observed trends and variability. Accordingly, continued efforts to improve model-data agreement in [O2] at the process level are required, in order to gain a mechanistic understanding of observed changes in oceanic oxygen towards improved predictions of future change. To this end, model experiments conducted here reveal that the sign and magnitude of [O2] change over the last approximately 50 years depends critically on ocean acidification feedbacks and prescribed air-sea fluxes of heat, water and momentum.

    Our results suggest that explicit representation of observationally based ocean acidification impacts on photosynthetic carbon drawdown yields major changes in the spatio-temporal distribution of simulated POC export and, consequently, subsurface O2 utilization. Historical POC export changes associated with the inclusion of a pCO2-sensitive C : N ratio are similar in magnitude but differ in sign to the absolute historical changes predicted by current fixed stoichiometry models as a result of secular ocean warming. For example, Laufkötter et al. [87] report POC export at 100 m decreases by approximately 0.5 mol C m−2 yr−1 for much of the extrapolar ocean between 1960 and 2006 in hindcast simulations of the CCSM3-BEC model. With the inclusion of a pCO2-sensitive C : N ratio, POC export at 100 m increases by more than 0.3 mol C m−2 yr−1, suggesting that stoichiometric effects could act to compensate a significant component of the ‘direct’ climate-driven reduction in export production. Thus, as highlighted by Tagliabue et al. [33], the inclusion of stoichiometric plasticity in the next generation of ESMs could alter the classical view that historical [87] and future [88] anthropogenic forcing drives a simulated reduction in global marine production, the observational basis for which remains contested [89,90]. Equally, Laufkötter et al. [87] report increases in POC export of up to 0.6 mol C m−2 yr−1 in subpolar (light-limited) regimes, associated with historical density stratification increases in their hindcast model. In this case, the inclusion of ocean acidification effects on organic carbon drawdown could act to amplify secular increases in POC export change associated with ocean warming. The development of high complexity ocean biogeochemistry models, which dynamically resolve stoichiometry (e.g. ERSEM [91]), provide a framework for investigating diverse ecosystem responses to changes in environmental conditions.

    A strengthened ‘soft tissue’ pump associated with the inclusion of a pCO2-sensitive C : N ratio also drives major changes in simulated remineralization within the ocean interior, as evidenced by stronger gradients in [DIC] and [O2]. Coeval subsurface [O2] decreases and a reduction in the global inventory of O2 are consistent with model projections which include pCO2-sensitive C : N ratios in carbon fixation (e.g. [32,33]). However, counter to published model experiments, the spatial pattern of [O2] decrease associated with carbon overconsumption in this study follows the absolute fingerprint of (observed) ventilation-driven deoxygenation [92], with the most pronounced O2 depletion occurring at mid- to high latitudes, rather than focused within the tropical OMZs. This result does not, therefore, support the suggestion (e.g. [4]) that inclusion of a pCO2-sensitive C : N ratio in carbon fixation could invoke a net deoxygenation of the tropical thermocline for the historical period, in agreement with recent hindcast [9] and prognostic [27] ESM experiments which include variable stoichiometry. Rather, stoichiometric effects in response to historical ocean acidification are shown here to bring about elevated oxygen depletion in mid- to high-latitude regions of water renewal, thus providing a biogeochemical amplifier for (chiefly) physically driven simulated [O2] changes (e.g. [6]). Accordingly, the inclusion of these effects in the next generation of ESMs could act to reduce the discrepancy between observed [91] and more muted ESM-simulated historical deoxygenation trends. Additionally, consistent with prognostic model studies, the inclusion of a pCO2-sensitive C : N ratio here produces an increase in the modelled volume of low-O2 waters. Although, due to the smaller imposed pCO2 forcing, increases in suboxic volume simulated for the historical period (+2% for 2013) are smaller than that projected for 2100 in variable stoichiometry models (+36–50%, [32,33]). Further expansion of suboxic and hypoxic water bodies in response to anthropogenic forcing has important implications for the global marine nitrogen cycle, with a coeval increase in simulated denitrification rates (34%; 0–2000 m) associated with the elevated consumption of nitrate (Inline Formula) in microbial decomposition of organic material under suboxic conditions. A series of parallel PlankTOM10-NEMOv.3.1 sensitivity experiments carried out using a diagnostic N2O model [73] also indicate that pCO2-sensitive C : N ratios drive enhanced historical marine N2O production via promotion of high-yield low-O2 processes (e.g. denitrification and nitrifier-denitrification).

    As noted by Matear & McNeil [93], however, caution is required when extrapolating the results of mesocosm experiments carried out with one natural plankton assemblage [28] to the global scale for all phytoplankton taxa and biogeographic provinces. Therefore, while model parametrizations based around the results of [28] comprise an aggregation of data points across nine large mesocosm enclosures at varying pCO2, further experimental evidence from diverse marine systems under different environmental conditions is required to better constrain this potentially significant ocean carbon cycle feedback in models. To this end, Hutchins et al. [94] highlight that, despite good agreement within unialgal cultures towards elevated C : N ratios in response to acidification, reported stoichiometric changes within CO2 manipulation experiments carried out on natural plankton communities are much more variable. Differences between results are attributed in part to changing experimental practice, such as ‘batch-mode’ versus continuous culture incubation methods. However, a number of external biological factors are also invoked to explain the inconsistent response of natural assemblages to elevated pCO2, including differing zooplankton grazing rates or community composition between oceanographic regimes. For instance, Thingstad et al. [95] demonstrate for the Arctic pelagic ecosystem that elemental stoichiometry changes within phytoplankton biomass depend critically on the nature of growth-limiting factors within the heterotrophic community. Additional uncertainties are entrained into biogeochemical model experiments which include variable C : N ratios due to the lack of a process-level understanding for the excess carbon uptake reported by Riebesell et al. [28] in response to rising pCO2. For instance, as suggested by Arrigo [96], additional fixed DIC could be released as dissolved organic carbon and contribute to carbon export through the formation of transparent exopolymer particles [97] via increased aggregation and particle sinking.

    Our results also suggest that accounting for pCO2-driven perturbation to biogenic calcification rates has only a negligible impact on model simulated [O2] dynamics for the historical period. However, as highlighted by Ridgwell et al. [98], unresolved questions regarding the observed ‘form and sensitivity’ of ocean acidification impacts on calcification introduce uncertainties into the parametrizations adopted in global models such as PlankTOM10. Moreover, the impact of ocean acidification on the ‘PIC : POC rain ratio’ must be considered alongside other external influences on the biological pump, particularly changes in export production in response to ocean warming [6,88]. For instance, model simulated historical [87] and future [6,88] reductions in overall export production associated with increased density stratification may act to moderate the impact of any acidification-driven reductions in CaCO3 production on Inline Formula. Although, while this study suggests that pCO2-calcification effects may be less important for biogeochemical cycles over the historical period, these processes may remain relevant on centennial [38] and millennial [42] time scales.

    Finally, our results suggest that imposed atmospheric forcing plays a major role in modulating interannual variability in subsurface [O2]. Particularly, in agreement with forced ocean biogeochemistry model results of [86] for the North Atlantic, simulated interannual variability in [O2] is shown here to be primarily associated with heat fluxes (σθ) in extratropical (mainly subpolar) regions and wind stress (σwspd) in the tropics. Large uncertainties between prescribed meteorological datasets in these regions propagate into simulated thermocline [O2], consistent with a recent intercomparison of surface reanalysis data, which attributes elevated multiproduct inconsistency in the tropics and extratropics chiefly to wind stress and heat flux uncertainties, respectively [99]. The important role of tropical zonal wind stress in controlling variations in model-simulated low-O2 water bodies has been demonstrated in recent work [24,100], providing further motivation for the provision of appropriate wind forcing to ocean-only models investigating O2 dynamics.

    While all atmospheric forcing products used here provide surface fluxes which are a priori representative of observed changes in meteorological variables (aside from IPSL), data choices still place major constraints on simulated changes in ocean properties. Further work is required to better understand the biases in atmospheric forcing datasets, both in terms of comparing meteorological fields (e.g. [99]) and assessing how uncertainties in these prescribed forcings (and bulk formulae) impact upon the evolution of hindcast variables. Towards this objective, the Coordinated Ocean-ice Reference Experiments (COREs) project proposes a standard protocol for running hindcast ocean-ice models, emphasizing the need for models to be integrated using different atmospheric forcings in order to ‘assess implications on the ocean and sea ice climate of various atmospheric reanalysis or observational products' [101]. However, the majority of recent physical [102] and biogeochemical [103] multi-model hindcast studies remain focused on investigating the implications of a common atmospheric forcing for a range of ocean models. This approach assumes that all data products provide an equally appropriate representation of historical changes in observed air-sea fluxes. Rather, as argued here, a multifaceted approach is required in order to better evaluate surface meteorological data products, involving both multi-model intercomparison under a common atmospheric forcing (e.g. COREv2-IAF; [102]) and ensembles of different atmospheric forcing experiments using a common ocean model.

    In this study, we have presented an analysis of processes using a relatively coarse global ocean biogeochemistry model, which shares some of the shortcomings on the representation of [O2] as similar models noted in §1. To ensure the model shortcomings had a limited effect on the process-based results presented here, we initialized the model with observations and removed the ensuing drift by isolating processes using the REF reference simulation. This methodology is unlikely to influence the main conclusions of our analysis on the potential influence of specific processes and the large role of atmospheric forcing for variability, but it could have a larger influence on N2O fluxes which are more regulated by narrow ranges in low [O2]. Further work could examine the interactions between physical and biogeochemical factors using a higher-resolution model, particularly to better quantify processes that influence oceanic N2O emissions.

    Data accessibility

    Model data can be made available by contacting O.A. Model restart files are available at http://opendap.uea.ac.uk:8080/opendap/greenocean/Restart. Forcing and CMIP5 data used in this study are publically available.

    Authors' contributions

    Conception and design of the study by all authors. Model development and testing by E.B., C.L.Q. and O.A. Model runs, analysis and first manuscript draft by O.A. All authors contributed to the interpretation of results and drafting of the final manuscript.

    Competing interests

    We declare we have no competing interests.

    Funding

    This work was supported by the European Union's Horizon 2020 research and innovation programme under grant agreement no. 641816 (CRESCENDO).

    Acknowledgements

    We thank C. Enright for programming support. Model experiments and analyses presented in this paper were carried out on the High Performance Computing Cluster supported by the Research and Specialist Computing Support service at the University of East Anglia. We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modelling groups for producing and making available their model output. For CMIP, the US Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led the development of software infrastructure in partnership with the Global Organization for Earth System Science Portals.

    Footnotes

    One contribution of 11 to a discussion meeting issue ‘Ocean ventilation and deoxygenation in a warming world’.

    Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3816628.

    Published by the Royal Society. All rights reserved.

    References