Shrinking of the Arabian Sea oxygen minimum zone with climate change projected with a downscaled model

1 Introduction

A land-locked northern boundary and limited ventilation to the north, combined with intense surface seasonal blooms lead to the formation and maintenance of one of the thickest, open ocean Oxygen Minimum Zones (OMZs) in the Arabian Sea (AS; Figures 1A, F). As dissolved O₂ is crucial for sustenance of most marine animal life, the presence of low-oxygen mesopelagic zones can potentially impart stress on the marine life (Bopp et al., 2013), resulting in significant reduction in the marine biodiversity in OMZs as compared with well oxygenated regions (Stramma et al., 2012).

Figure 1 Evaluation of modelled present day O₂ against Word Ocean Atlas observations. The left column shows the O₂ inventory averaged over the layer 200-700m, the right column shows the meridional section at 65°E (indicated by black dashed lines on panel A) with white contours indicating the hypoxic regions (O2< 60 μmol/L). (A-F) WOA observations (WOA 2018); (B–G) downscalled experiment (ROMS_CTL); (C–E, H–J) the three CMIP5 models (MPI, GFDL, IPSL) over the present-day period (CMIP5_HIST simulations, averaged from years 1986 to 2005). The northern Arabian Sea (NAS) and southern Arabian Sea (SAS) regions are delimited by red dashed lines on panel (A).

The OMZs also play crucial role in the cycling of gases (e.g., CO₂, N₂O, CH₄), and elements (e.g., nitrogen, sulfur and carbon) that have direct implications for the Earth’s climate. For instance, the emission of potent greenhouse gases (such as nitrous oxide, methane, etc.) resulting from the anaerobic metabolism in the OMZs (Dileep, 2006) can exacerbate the greenhouse warming and the associated climate change (Naqvi & Noronha, 1991; Codispoti et al., 2001; Chronopoulou et al., 2017). In addition to this, the denitrification and anaerobic ammonium oxidation (anammox) in the OMZs lead to fixed nitrogen loss in the pelagic ocean (Stief et al., 2017), limiting the inventory of bioavailable oceanic nitrogen (Falkowski, 1997). This can potentially inhibit the marine productivity, and hence affect the marine food web. Another important concern arises from the direct interaction of open ocean OMZs with the neighboring coastal regions. For instance, the upwelling along west coast of India during boreal summer triggers the development of largest natural seasonal hypoxic (O₂< 60 μmol/L) events, causing frequent major fish-kill events there (Naqvi et al., 2009; Parvathi et al., 2017a; Pearson et al., 2022). On the other hand, the seasonal upwelling in this region bears a large sector of the Indian economy by supporting the pelagic fisheries during the summer monsoons (Vivekanandan et al., 2005). The recurring seasonal hypoxia and their sporadic intensification thus adversely affects the ecosystem and fisheries stock in this region and forms a direct threat to the Indian economy (Naqvi et al., 2009).

Despite the strong seasonal modulation of both atmospheric forcing and oceanic productivity, the AS OMZ is reported to have no pronounced seasonal changes (Resplandy et al., 2012). In the present-day climate, this has reportedly been maintained through the subtle balance that exists between supply of oxygen through dynamical processes (i.e., transport through western boundary currents, mixing and advection by mesoscale eddies and convective mixing), and the biological demand resulting from seasonal blooms (Sarma, 2002; Resplandy et al., 2012; McCreary et al., 2013). However, a recent modeling study by Lachkar et al. (2018) have demonstrated that both the source (dynamical supply) and sink (biological demand) of oxygen in the AS OMZ are highly sensitive to the monsoonal winds. In addition to this, the future projections from CMIP5 simulations indicate significant changes in both the dynamical forcing (atmospheric circulation) and associated biological productivity in response to climate change (Sandeep & Ajayamohan, 2015; Praveen et al., 2016; Roxy et al., 2016; Parvathi et al., 2017b). Moreover, the steady and heterogenous warming of AS reported at decadal to multi-decadal time-scales (Han et al., 2014; Roxy et al., 2016; Gopika et al., 2020; Lachkar et al., 2021) has the potential to reduce the solubility of O_2, and also strengthen the oceanic stratification, thus weakening the vertical mixing and limiting the supply of O₂ to the sub-surface layers. Additionally, the shallow marginal seas in the neighborhood of AS play a crucial role in controlling the extent of sub-surface OMZ (e.g., Naqvi et al., 2009; McCreary et al., 2013; Lachkar et al., 2019). The ongoing warming increases the buoyancy of the shallow Arabian Gulf (also known as Persian Gulf) water, preventing its subduction into the northern AS, and weakening one of the ventilation pathways of the AS OMZ leading to an expansion of its core (Lachkar et al., 2019; Lachkar et al., 2021). Despite all of the above facts, there is no clear consensus on the impact of climate change on the AS OMZ. The complex interaction between physical and biological processes and their differential responses to the climate sensitivity pose a major constraint for the accurate anticipation of projected response in the OMZ under climate change. However, most of the earlier findings indicate that despite its occurrence throughout the geological time-period (Fennel & Testa, 2019), the recent intensification of the global OMZs (both coastal and open ocean) is vastly due to the dominance of anthropogenic climate change (Breitburg et al., 2018). The important role of oxygen in maintaining a healthy ecosystem, diverse fisheries, and hence the economy of the highly populated coastal belt along the western India calls for a thorough assessment of the impact of climate change on the AS OMZ.

Climate simulations from the Earth System Models (ESMs) have widely been utilized to study the impact of climate change on a broad range of physical-biogeochemical systems as well as for making policies that help to mitigate the adverse impacts. While the ESM simulations performed under CMIP5 bring great opportunity to explore future projections of both atmospheric and oceanic processes, they also possess biases that can potentially alias the interpretation of climate change driven future projections (Lachkar et al, 2021; Gopika et al., 2020). Moreover, the biases are generally reported to be larger (in comparison to the physical parameters for instance) for biogeochemical variables, and particularly for oxygen (Bopp et al., 2013; Cocco et al., 2013; Schmidt et al., 2021). Earlier studies focusing on the north Indian Ocean/AS have indicated the cascading effects of misrepresentation of key processes in bringing huge biases in the mean-state of both atmospheric and oceanic simulated fields. For example, the misrepresentation of main water masses leads to the biases in depth of mixed layer/thermocline in the northern AS, which in turn affects the present-day simulations of Indian monsoons in the CMIP5 models (Nagura et al., 2018). The misrepresentation of water masses in the CMIP5 models have also been identified as the main reason for the biases in the AS OMZ simulations by Schmidt et al. (2021). The ill-represented AS OMZ mean state in the CMIP5 prevents utilizing the ensemble of these simulations alone to draw meaningful conclusions in a regional context (Schmidt et al., 2021).

Another important contribution for the supply of oxygen to the interior of the ocean and hence the ventilation of open ocean OMZs comes from the mixing (both vertical and horizontal transport; Lévy et al., 2021). The OMZ in the AS is largely influenced by the strong mixing in its western boundary arising from energetic mesoscale eddies (McCreary et al., 2013). Hence, the simulation of present-day mean state distribution of oxygen in the AS demands an accurate representation of the mesoscale spectrum, either through increasing grid resolution or though choosing accurate parameterization schemes (Resplandy et al., 2012; Lachkar et al., 2016). The computational cost and the global coverage prevent the ESM under CMIP5 to account for such great details in the regional basins. In order to overcome some of the biases of global CMIP5 simulations possibly arising from the above factors, and gain some insight about the robust patterns of projected changes in AS OMZ in response to climate change, along with its driving processes, the present study, attempts a regional downscaling approach. Such downscaling methods have been used to understand the physical-biogeochemical interaction in the future climate change scenario elsewhere in the global oceans (Echevin et al., 2020; Buil et al., 2021), but not yet in the Indian Ocean.

Hence in the present study, we specifically address the following questions: (1) Will the AS OMZ expand or shrink in response to climate change? (2) What are the contributions of local vs. remote forcing in driving the projected O₂ changes? (3) How sensitive are future response patterns to the model representation of the present-day state? To this end, we performed a set of simulations with an eddy-permitting (1/3°) regional configuration of the physical-biogeochemical model ROMS over the AS domain, using atmospheric fields and oceanic boundary conditions taken from observed climatologies and from CMIP5 model projections. We have designed a set of specific sensitivity experiments to identify the evolution of key physical and biogeochemical processes that maintain the present-day balance of the AS OMZ. Finally, we show the potential and the limits of regional downscaling to enhance our ability to predict future OMZ changes and more generally improve our understanding of the key biogeochemical processes in the regional basins.

The paper is organized as follows: Section 2 is dedicated to provide the details of the methodology. In particular, this section includes the details of the CMIP5 simulations, regional model experiments with ROMS, and observations. Section 3 describes the main results of our study. In section 4, we provide a discussion of the lessons learnt and the limitations of the study. Finally, section 5 provides the main conclusions of the study.

2 Methods

In this section, we introduce the three selected CMIP5 models and the associated downscaled eddy-permitting ROMS experiments, followed by details of our downscaling approach. A summary of these details is provided in Tables 1, 2.

Table 1 Description of the three CMIP5 models used in the present study.

Table 2 Description of the regional, downscalled simulations of the Arabian Sea performed with the ROMS model, and of their forcings.

2.1 CMIP5 simulations

The CMIP5 archive contains global coupled simulations performed with Earth System Models (ESM). These simulations cover the present-day “historical” period (denoted as CMIP5 _HIST_modelname in this study), where the coupled model is forced by observed atmospheric composition changes of greenhouse gases during the industrial period (mid-nineteenth century to near present). They are followed by future projection simulations, known as Representative Concentration Pathway (RCP), which are forced with atmospheric composition changes consistent with varying degrees of greenhouse gas emission scenarios (denoted as CMIP5 _RCP_modelname in this study). The RCP experiments are named based on the amount of radiative forcing added to the experiment. Here, we used the unmitigated high emission scenario RCP8.5 wherein the radiative forcing increases throughout the 21^st century before reaching a level of 8.5 W m^-2 at the end of the century (Taylor et al., 2012; Bopp et al., 2013).

We focus our analysis on a subset of simulations from the CMIP5 archive, based on three ESM, i.e., MPI-ESM-LR, GFDL-ESM2M, and IPSL-CM5A-MR (Table 1, here onwards we refer them as MPI, GFDL and IPSL). These three ESM are based on different ocean general circulation models, along with different mixing schemes and resolutions (both horizontal and vertical), and on different biogeochemical models, with various degrees of complexity. They also differ in representing the influence of marginal seas, particularly of the Arabian Gulf (AG). For instance, the MPI model includes the AG water explicitly, and uses the Pacanowski and Philander (1981) mixing scheme. In the GFDL model, the impact of AG water is parameterized by enhanced cross-land mixing. The GFDL model uses the KPP mixing scheme. The IPSL model does not include the effect of AG water explicitly nor does it represent its impact through restoring temperature or salinity to observation. The mixing scheme used in the IPSL model is Blanke and Delecluse (1993) (A detailed discussion on this can be found in Huang et al., 2014 and Nagura et al., 2018).

These differences explain in part the different mean states of the OMZ represented in the historical simulations (Figure 1). In fact, this subset of three CMIP5 models was actually selected (out of 9 available) because it covers contrasted representations of the present-day OMZ (Figure 1). CMIP5_HIST_MPI (Figures 1C, H) has the closest agreement with observations (Figures 1A, F), CMIP5_HIST_GFDL (Figures 1D, I) has intermediate skill and CMIP5_HIST_IPSL (Figures 1E, J) completely fails at capturing the observed subsurface depletion of oxygen in the AS, however has lesser O₂ in the sub-surface layers compared to the surface (Figure 1J). The representation of the OMZ in each simulation is quantified more accurately by comparing the volume fraction of the OMZ (defined as the volume of hypoxic waters, i.e. where O2<60 μmol/L in the region north of Equator and comprised between 0 and 1000 m depth) and of the core of the OMZ (defined as the volume of suboxic waters, i.e. where O2<4 μmol/L in the same region) in the simulations with climatological observations (Figure 2). Such quantitative assessments are important to determine the reliability of model’s skill in representing the observed variables. Our analysis reveals that, CMIP5_HIST_MPI, which best captures the volume of the OMZ, does it at the expense of a strong overestimation of the volume of suboxic waters (which was not very clear from the qualitative maps), while the CMIP5_HIST_GFDL does not have any suboxic waters, and the CMIP5 _HIST_IPSL does not even represent hypoxic waters (Figure 2).

Figure 2 Quantitative evaluation of modelled present day OMZ against Word Ocean Atlas observations. Volume fraction occupied by suboxic water masses (OMZ_core, O2< 4 μmol/L) and by hypoxic water masses (OMZ, O2<60 μmol/L) in the 0 – 1000 m top layer of the Arabian Sea, in Word Ocean Atlas observations (WOA 2018, black), in CMIP5 historical simulations (shades of red: CMIP5_HIST_MPI, CMIP5_HIST_GFDL, CMIP5_HIST_IPSL), and in the downscalled ROMS present-day simulation (ROMS_CTL, in blue). Note that the only CMIP5 model that simulates OMZ_core waters is the MPI model, and that the IPSL CMIP5 model does not even simulate OMZ waters.

2.2 Regional simulations

2.2.1 Model configuration

We used the Regional Ocean Modeling System Adaptive Grid Refinement in Fortran (ROMS-AGRIF; http://www.croco-ocean.org/; here onwards mentioned as ROMS) with an embedded biogeochemical model to simulate the dynamics of the AS-OMZ (Lachkar et al., 2016). ROMS has terrain following sigma coordinates and solves the primitive equations with Boussinesq and Hydrostatic approximations. (Shchepetkin and McWilliams, 2005). To reduce spurious mixing, ROMS includes a rotated-split third-order upstream biased algorithm for the advection of tracers (Marchesiello et al., 2009). The subgrid vertical mixing is represented by including non-local K-profile parameterization (KPP) scheme (Large et al., 1994).

The embedded biogeochemical model has previously been used in ROMS simulations of the AS OMZ (Lachkar et al., 2016; Lachkar et al., 2018; Lachkar et al., 2021) and is fully described therein. It includes nitrate, ammonium, one class of phytoplankton, one class of zooplankton, two classes of detritus, a dynamic chlorophyll-to-carbon ratio and a module that describes the evolution of O₂, with anaerobic respiration below a fixed threshold (O₂< 4μmol/L) and benthic denitrification.

Our regional configuration covers the AS meridionally from 5° S to 30° N and zonally from 33° E to 78° E, with lateral open boundaries at 5° S and 78°E. The model is forced at the surface with wind, air-temperature, humidity, SST and atmospheric heat fluxes, and with temperature, salinity and currents at lateral boundaries. Apart from the above, the oxygen and nitrate initial and boundary conditions were also provided for the biogeochemistry. Depending on the experiment, the surface forcing and lateral boundary fields were either constructed from observations or from CMIP5 outputs (details are provided in the subsequent sections). Experiments were performed using an eddy-permitting 1/3° horizontal grid resolution with 32 sigma layers. This allowed a significant improvement in the representation of the Arabian Sea OMZ relative to the CMIP5 models (Figures 1, 2). Using the downscaled model at higher horizontal resolutions (1/6° and 1/12°) brought about less changes in terms of the simulated present day OMZ than changes associated with the model (CMIP5 vs ROMS, Figure SI 1). Therefore, for computational efficiency and to be able to run multiple sensitivity simulations, we decided to use this same resolution (1/3°) for all regional downscaling experiments presented in this study. The sensitivity to model resolution of the response of the OMZ to future climate change will be the topic of a follow-up study. Contrary to CMIP5 projections, we have the same base model for all the downscaled experiments, which is one of the main advantages of our downscaling approach. Hence, this approach will allow us to examine the susceptibility of future projections to model’s representation of the present-day state and the contributions of local and remote forcing to projected changes.

2.2.2 Control present-day experiment (ROMS_CTL)

The control simulation (denoted as ROMS_CTL throughout the paper) is forced with observation-based climatologies, and thus represents the present-day (equivalent to CMIP5 historical) period. The atmospheric forcing fields for heat and freshwater fluxes are constructed from COADS monthly climatology, and we used monthly wind stress climatology based on QuikSCAT - derived Scatterometer Climatology of Ocean Winds (SCOW; Risien and Chelton, 2008). Temperature, salinity, oxygen, and nitrate were initialized and forced at the open boundaries using monthly climatology from World Ocean Atlas (2018). Currents were initialized and laterally forced at the open boundaries from Simple Ocean Data Assimilation (SODA) at monthly frequency. For all our simulations, we applied a restoring of surface temperature and salinity using kinematic flux corrections described in Barnier et al. (1995). The ROMS_CTL experiments were run for 50 years. The model was spun-up until a steady state was reached (45 years) and the last 5 years were used for analysis (Figure SI 2).

The overall evaluation of physical and biogeochemical properties in summer and winter in the ROMS_CTL experiment is shown in Figures SI 3-5. The main features of the model solution are the signature of seasonal upwelling along the coasts of Somalia, Oman and western India comparable with observations, with however a slight overestimation of surface chlorophyll in the Great Whirl region. The winter convective mixing in the northern AS is also reasonably captured as can be inferred from the cooling of SST (Figure SI 3), along with nutrient injection (Figure SI 4) and enhanced surface chlorophyll (Figure SI 5). Our model also captures the main circulation patterns and its seasonal modulations in the AS domain (not shown). However, from the water mass analysis shown in Figure SI 6A, E), we also notice that, the ROMS_CTL do not simulate the high salinity tongue - indicating the underestimation of Gulf water intrusion in the northern AS (our model domain includes the Arabian Gulf). This underestimation was already reported in Lachkar et al. (2018), and could be due to the model limitations arising from the lack of realistic forcing data in the Gulf region or the model grid resolution (Lachkar et al., 2018).

Despite the above factors, the comparison of the volume fractions of OMZ and OMZ_core in the AS from the ROMS_CTL experiment with observation reveals better skills of the ROMS_CTL experiment at capturing the AS OMZ, both its core and total volume, compared to the set of CMIP5_HIST simulations (Figure 2).

2.2.3 Regional downscaling experiments

2.2.3.1 Direct downscaling (DD)

In direct downscaling (DD) approach, we use CMIP5 model simulations to initialize and force our regional AS simulations with no correction of the CMIP5 model biases (with respect to observations). We applied DD to both historical (denoted as ROMS_HIST) and future (referred to as ROMS_DD) simulations. For historical simulations, we used monthly climatologies from historical CMIP5 simulations (i.e., CMIP5_HIST) computed over 1986-2005 as forcing and boundary (for physical variables as well as biogeochemical tracers) fields. Similarly, for future simulations, we used monthly climatologies from future CMIP5 simulations (i.e., CMIP5_RCP8.5) averaged over 2081 to 2100. As for ROMS_CTL, both historical and future climate downscaled experiments were run for total 50 years with a spin-up-up of 45 years. We used the last 5 years for the analysis presented here.

2.2.3.2 Bias corrected downscaling (BCD)

An additional downscaling approach, denoted as Bias-Corrected Downscaling (BCD) and similar to the Pseudo Global Warming downscaling method discussed in earlier regional climate studies (Xu et al., 2019), is used here. In this approach we correct the present-day biases in CMIP5 models with respect to observations before using this data to force the future regional projections (referred to as ROMS_BCD). More concretely, we added future anomalies from CMIP5 models to observational climatologies (used to force ROMS_CTL and denoted OBS) to construct the forcing and boundary fields (i.e., (CMIP5_RCP - CMIP5_HIST) + OBS) used to initialize and force the future regional AS simulations (i.e., ROMS_BCD). The monthly averages of CMIP5 future anomalies were computed by subtracting the last 20 years averages of ESM simulated variables for present-day (i.e., 1986-2005, CMIP5_HIST) from the future simulations (i.e., 2081-2100, CMIP5_RCP; Figure SI 7, 8, only O₂, and NO₃ are shown). As for the experiments mentioned in previous sections, these downscaling experiments were also run for 50 years, where the model was spun-up for 45 years and, the last 5 years were used for the analysis. In both DD and BCD, when the monthly outputs were not available (for O₂ & NO₃), we repeated the annual values to construct the monthly cycle.

2.2.3.3 Quantification of the contributions of local vs remote forcing to Ochanges:

To disentangle the roles of local vs. remote forcing in controlling the future evolution of the AS OMZ, we performed a first set of sensitivity downscaling experiments (referred to ROMS_BCD_LOC) where the prescribed atmospheric forcing is similarly based on bias-corrected future conditions (similar to ROMS_BCD runs) while lateral boundary conditions are maintained at the current conditions (similar to ROMS_CTL). Additionally, to explore the role of physical vs. biogeochemical forcing we performed a second set of twin simulations (denoted ROMS_BCD_PHY), where only biogeochemical boundary conditions (i.e., O₂, NO₃) are maintained at current conditions (similar to ROMS_CTL) whereas both atmospheric and physical boundary conditions (i.e., T, S, u, v) are based on future conditions (similar to ROMS_BCD runs).

Using our sensitivity simulations, we linearly decompose the downscaled projected O₂ changes obtained by subtracting the control simulation ROMS_CTL from the bias-corrected downscaled experiment (Total changes = ROMS_BCD – ROMS_CTL) into the sum of three contributions: (i) local forcing, (ii) remote physical forcing and (iii) remote biogeochemical forcing.

While the effect of changes in local forcing on future O₂ changes (${\text{LOC}}_{\text{FRC}}$) can be assessed by simply subtracting the ROMS_CTL from the bias-corrected downscaled experiment forced with present-day boundary conditions ($\text{i}.\text{e}.,\text{ ROMS}\_{\text{BCD}}_{\text{LOC}}-\text{ROMS}\_\text{CTL}$), the impact of changes in physical conditions at the lateral boundaries can be evaluated by subtracting the BCD experiments forced with present-day boundary conditions ($\text{ROMS}\_{\text{BCD}}_{\text{LOC}})$ from the BCD experiments forced with present-day biogeochemical boundary conditions ($\text{ROMS}\_{\text{BCD}}_{\text{PHY}}$ ). Finally, the contribution of changes in biogeochemical conditions at the lateral boundaries (${\text{REM}}_{\text{BGC}}$ ) can be assessed by subtracting the BCD experiments forced with present-day biogeochemical boundary conditions from the BCD experiments forced with future conditions ($\text{i}.\text{e}.,\text{ ROMS}\_\text{BCD}-\text{ROMS}\_{\text{BCD}}_{\text{PHY}}$).

Following earlier studies (e.g., Lachkar et al., 2020), we assessed the impact of thermal and non-thermal drivers on the total O₂ changes from both CMIP5 and downscaled simulations. For this, we decomposed the total O₂ changes into O₂ saturation (O₂Sat), accounting the thermal effects, and Apparent Oxygen Utilization (AOU), accounting the non-thermal effects (i.e., ventilation as well as biological components, Kwiatkowski et al., 2020).

3 Results

3.1 Future O₂ changes in the AS OMZ

We analyze the impact of climate change on the AS OMZ (O₂ averaged over 200-700m layer) by contrasting the future projections from both CMIP5 and downscaled simulations (Figure 3). All CMIP5 models consistently project future oxygenation (positive anomaly, Figures 3A–C) of the sub-surface layer in the southern AS (SAS, Equator to ~10°N) in agreement with earlier studies (e.g., Bopp et al., 2013; Cocco et al., 2013; Bianchi et al., 2018). A similar oxygenation projected from all 3 models in two regional downscaled approaches further strengthens the reliability of this O₂ future response pattern (Figures 3A.1-C.1, A.2–C.2). In contrast to this, both the magnitude and sign of projected O₂ in the northern AS (NAS, north of ~20°N) exhibits a strong sensitivity to the choice of CMIP5 model (Figures 3A–C). The deoxygenation is stronger and encompasses larger area over the NAS in CMIP5_GFDL (Figure 3B) compared to much weaker deoxygenation seen in both CMIP5_MPI (Figure 3A) and CMIP5_IPSL (Figure 3C). These differences may in part reflect differences in the representation of local physics and biogeochemistry in these models. To explore this hypothesis, we contrast regional downscaled models with their parent CMIP5 model. Using a common physical-biogeochemical model (ROMS-NPZD), the DD projections expectedly display slightly weaker model discrepancies (associated with differences in parent model present-day conditions) relative to those in CMIP5 models. This is particularly true in the NAS where the strongly contrasting O₂ changes simulated by the three CMIP5 models are significantly reduced in the DD experiments (Figure 3). However, despite using a common physical and biogeochemical model, these simulations still display substantial differences in their projected O₂ changes in both NAS and SAS. This indicates that, despite the higher regional model resolution and the likely improved representation of local physics and biogeochemistry in the regional simulations, other factors and biases inherent to the parent models such as the representation of present-day conditions and their projected large-scale physical and environmental changes strongly influence the downscaled O₂ projections.

Figure 3 Projected oxygen changes in the OMZ layer (200-700m). (A–C) with the 3 CMIP5 experiments; (A.1–C.1) with the corresponding Direct Downscalling ROMS experiments; (A.2–C.2) with the corresponding Bias-Corrected Downscalling ROMS experiments. Oxygen losses (deoxygenation) appear in blue shading, and oxygen gain (oxygenation) in red shading. Red dots show the location of anoxic waters (OMZ core) and black dots show the location of hypoxic waters (OMZ) in the corresponding present-day simulations (CMIP5_HIST for CMIP5, ROMS_HIST for DD, ROMS_CTL for BCD). By construction, the present-day experiment (ROMS_CTL) is the same for the three BCD experiments, and the BCD experiments are the best estimates.

3.2 Sensitivity of future oxygen projections to model representation of present-day conditions

The DD and the BCD simulations differ only by their initial state, as they are conducted with a common (ROMS-NPZD) model, and are forced with the same CMIP5 anomalies at the boundaries. Therefore, contrasting the projections computed from these experiments allows us to quantify the sensitivity of the future oxygen projections to model representation of present-day conditions. The comparison between the BCD and DD projections indeed reveals significant differences in the OMZ layer in NAS, particularly with models that are associated with large biases in their representation of the present-day OMZ, such as the GFDL and IPSL (Figure 3). However, when corrected for the present-day biases, the spread in the projected patterns in the NAS reduces substantially among the three models in BCD experiments relative to the DD simulations (Figure 3). Indeed, the three BCD simulations consistently project a slight oxygenation in the NAS in opposition with the conflicting trends in the DD and CMIP5 simulations. In contrast to this, the inter-model differences are maintained and remain similar between the DD and BCD experiments in the SAS (Figure 3). These differences can also be seen in the vertical profiles of the projected O₂ anomalies presented in Figure 4. By removing the present-day biases in their initial state, the 3 models agree on the simulated changes throughout the water column in the NAS (Figure 4C). Moreover, in contrast to both CMIP5 (Figure 4A) and DD (Figure 4B) simulations, the layer of maximum deoxygenation is much shallower in the BCD experiments in the NAS (Figure 4C). In the SAS, however, the model discrepancies rather than the experiment strategies remain important among the BCD experiments and are similar to those simulated in the CMIP5 and DD runs (Figure 4). This indicates that the sensitivity of the O₂ projections to model representation of present-day conditions is much more important in the NAS relative to the SAS. It also suggests that downscaling following the BCD approach can help reduce model discrepancies and narrow down uncertainties in projected future O₂ changes in the NAS, while it adds little value in the SAS. Our quantitative analysis (Figure 5) further reveals that, the consistent oxygenation of the sub-surface layers in SAS leads to a future shrinking of AS OMZ under the climate change.

Figure 4 Projected changes in O2 vertical profiles. O2 response profiles averaged over the northern (NAS) and southern (SAS) Arabian Sea (A, D) with the three CMIP5 experiments, (B, E) with the corresponding Direct Downscaling ROMS experiments and (C, F) with the corresponding Bias-Corrected Downscaling ROMS experiments (MPI in blue, GFDL in red and IPSL in green). The NAS and SAS zones are shown in Figure 1A.

Figure 5 Projected changes in OMZ volume. Projected changes in the Volume fraction occupied by suboxic water masses (OMZ_core, O2 < 4 μmol/l) and by hypoxic water masses (OMZ, O2 <60 μmol/l) in the 0 – 1000 m top layer of the Arabian Sea, (A) in the three CMIP5 experiments, (B) in the corresponding Direct Downscaling experiments and (C) in the corresponding Bias-Corrected Downscaling experiments (MPI in blue, GFDL in red and IPSL in green). The OMZ volume is shrinking in all experiments.

In the next section, we examine the role of local vs. remote forcing in causing O₂ changes in both NAS and SAS regions.

3.3 Contribution of local vsremote forcing to projected future Ochanges

As both local and remote forcings play major role in controlling the coastal and open-ocean biogeochemistry in AS (Currie et al., 2013; Parvathi et al., 2017a; Pearson et al., 2022), we further use our downscaled model to disentangle their respective contributions to the future projected OMZ response. In Figure 6, we decompose the total projected O₂ response in the OMZ layer in to that resulting from local atmospheric changes, and remote physical and biogeochemical changes. This decomposition reveals that, each process contributes significantly to the total O₂ changes in the AS OMZ layer, depending on the region and the model (Figure 6). In the SAS, local forcing is leading to a weak deoxygenation along the Indian west coast and oxygenation over the rest, consistently among all models (Figures 6A.1-C.1). In all the models, the oxygenation along the west coast of India (region within SAS) is driven by both remote physical and biogeochemical processes, while it is dominated by the remote biogeochemical forcing in the rest of SAS (Figures 6A.2–C.2, A.3–C.3).

Figure 6 Projected oxygen changes in the OMZ layer (200-700m) (A–C) The impact of climate change on the ASOMZ computed from the 3 Bias-Corrected Downscalled ROMS experiments (ROMS_BCD – ROMS_CTL), (A.1–C.1) with a set of sensitivity experiments with which the impact of changes in surface forcing is isolated (ROMS_BCD_LOC − ROMS_CTL), (A.2–C.2) with a set of sensitivity experiments with which the role of changes in physical lateral boundary condition is isolated (ROMS_BCD_PHY − ROMS_BCD_LOC ), (A.3–C.3) with a set of sensitivity experiments the impact of changes in biogeochemical lateral boundary conditions are isolated (ROMS_BCD − ROMS_BCD_PHY). Oxygen losses (deoxygenation) appear in blue shading, and oxygen gain (oxygenation) in red shading. Differences between experiments reveal the respective contribution of local and remote forcings to the total projected change in O2.

3.4 Mechanisms driving the future OMZ changes

In the subsequent analysis, we separated out the contribution from thermal and non-thermal components on the O₂ changes (Figure 7). The nearly compensating influence of thermal and non-thermal components in the sub-surface layer of NAS leads to weaker OMZ changes in all 3 BCD projections (Figure 7). In contrast, in the SAS the non-thermal component, i.e., AOU, drives most of the sub-surface oxygenation (Figure 7). As the AOU is affected by both ventilation and biological changes (Kwiatkowski et al., 2020), the oxygenation (deoxygenation) signal can be indicative of either an increase (decrease) in dynamical supply or decline (increase) in biological demand. Within the scope of current study, we test the dominant driver for OMZ changes by analyzing the future projected annual productivity simulations from both CMIP5 and downscaling experiments (Figure 8). In agreement with earlier studies (e.g., Praveen et al., 2016; Roxy et al., 2016), all the models from CMP5 and regional downscaling project a consistent weakening of annual productivity in large parts of the AS under the climate change. The consistent weakening of biological demand arising from this decline in productivity in SAS (Figure 8) along with the larger supply of O₂ (Figure SI 7) at the southern boundaries, lead to a future oxygenation of the SAS under climate change. Whereas, in the NAS, the projected annual productivity response from both CMIP5 and regional downscaling is highly sensitive to the model.

Figure 7 Contribution from thermal and non-thermal components to the projected OMZ changes. (A-C) Projected changes in O2 averaged over the 200-700m layer in the Arabian Sea in the 3 Bias-Corrected Downscaling ROMS experiments, decomposed into their thermal (A.1-C.1; ΔO2sat) and non-thermal (A.2-C.2; - ΔAOU) components.

Figure 8 Drivers of oxygenation at the SAS. (A–C) Projected changes in the annual net primary productivity from 3 CMIP5 experiments; (A.1–C.1) with the corresponding Direct Downscalling ROMS experiments; (A.2–C.2) with the corresponding Bias-Corrected Downscalling ROMS. A future weakening appears in blue shading, and strengthening in red shading.

4 Discussion

4.1 Lessons learnt from the regional downscaling:

(a) Will the AS OMZ expand or shrink under the climate change?

The subtle balance between the dynamical ventilation and biological demand is reported to maintain the AS OMZ with no pronounced seasonality in the present-day climate (Sarma, Sarma, 2002; Resplandy et al., 2012 and the references therein). Previous studies based on CMIP5 simulations have exhibited the potential impact of climate change on both atmospheric circulation and the oceanic primary productivity (Sandeep & Ajayamohan, 2015; Praveen et al., 2016; Roxy et al., 2016; Parvathi et al., 2017b). Those projected changes have the potential to impact the AS OMZ in future. Using regional downscaling of a subset of ESMs under CMIP5, we find a consistent oxygenation in the sub-surface (200-700m) layer of the SAS, that potentially leads to shrinking of the AS OMZ (Figure 5) under climate change. This is strikingly different compared to other low-oxygen zones of the global oceans, which are reported to expand under the climate change (e.g., Keeling et al., 2010).

(b)How does the model representation of present-day conditions influence the future projections?

The convective mixing driven by northeast (boreal winter) monsoon winds acts as a ventilation pathway for the perennial AS OMZ (Morrison et al., 1999). This density driven overturning is highly influenced by the influx of high saline water from the Arabian Gulf. Moreover, the Arabian Gulf water itself acts as a ventilation pathway for the AS OMZ (Naqvi et al., 2009; Lachkar et al., 2019; Lachkar et al., 2021).

Our analyses indicate that misrepresentation of the main water masses and mixing in the northern AS in CMIP5 models has led to intense deoxygenation in those biased models. For instance, overestimation of inflow of the high saline Gulf water and the Arabian Sea High Salinity Water results in a very weak stratification in both GFDL and IPSL models (Figure 9). This weak stratification leads to anomalous mixing and hence deepening of the mixed layer, allowing intrusion of high oxygenated water into greater depths in the northern AS (Figure 1). However, the basin-wide warming (~4°C, calculated from 3 CMIP5 models) and the reduced winter monsoonal circulation (Parvathi et al., 2017a) under climate change act to weaken the vertical mixing, leading to intense deoxygenation in the future projections of GFDL. The above present-day misrepresentations have led to inconsistent future projections in the NAS. On the other hand, we have consistent response in the NAS when corrected for the present-day biases (Figures 3, 4). When the background stratification is very strong (Figure 9, the pattern is very similar even with salinity or density), even if we force the regional (ROMS_DD GFDL) model with biased surface and boundary fields, the projected deoxygenation in the north remains much weaker compared to the its parent CMIP5 model (Figure 3). This clearly indicates the sensitivity of mean state of model in determining the future projections of OMZ in the northern AS. When corrected for the present-day biases in regional models, we find that the layer of deoxygenation has shown significant shift and agreement amongst the models (Figure 4). This is particularly important as the deoxygenation and the associated shoaling of low oxygen water are known to impart stress on the ecosystem leading to habitat compression – a major negative impact of climate change as reported in many earlier studies (Eg., Köhn et al., 2022).

Figure 9 Role of background stratification in determining the vertical distribution. Vertical section of annual temperature (°C; in color shading), dissolved oxygen (μmol/L; dashed contours from 80 μmol/L and above concentration) along with winter mixed layer (m; asterisks) averaged in the northern Arabian Sea from observation (upper top, WOA2018), (A–C) the three CMIP5 models (MPI, GFDL, IPSL) over the present-day period (CMIP5_HIST simulations, averaged from years 1986 to 2005); (A.1–C.1) the corresponding (to CMIP5_HIST) Direct Downscalling ROMS experiments; and (A.2–C.2) from the corresponding (to CMIP5_HIST) Bias-Corrected Downscalling ROMS experiments.

(c) What drives the robust oxygenation in the southern Arabian Sea?

Earlier studies have shown the importance of southern boundary in ventilating the AS OMZ and thereby preventing its extension to the south (Sarma, 2002; Resplandy et al., 2012; McCreary et al., 2013). Using specific sensitivity experiments, we investigated the importance of southern boundary for the future evolution of AS OMZ under climate change. Our analysis revealed that, the oxygenation in the sub-surface layers of southern AS is largely driven by the remote processes influencing the southern and eastern boundaries. The reduced demand resulting from weaker primary productivity overwhelms the weaker dynamical supply under the climate change. With specific sensitivity experiments, we also demonstrate that the transport of O₂ from the southern Indian Ocean is important even under the climate change.

(d) How important is it to do regional downscaling for future projection of OMZs?

Unlike many parts of the global ocean, climate change induces an oxygenation of the southern boundary of AS that can potentially lead to a future shrinking of AS OMZ. From our downscaling experiments, we find that regional downscaling of global models is relevant only where local forcing is dominant, away from boundaries. For instance, in regions like northern AS, where there are many controlling factors for the present-day state, it is important to account for the mean-state biases. A correction of such biases can bring completely new insights, and also improve the inter-model agreement. The subtle balance between the source and sink can be highly sensitive to the layer and the region which we consider. With present study, we also learnt that, the commonly accepted practice of considering multi-model mean from CMIP5 can potentially mislead, as some of the models with anomalous response pattern can completely alias the mean response patterns. Hence, one has to practice caution while interpreting the impact of climate change on the dissolved O₂ from global coupled models.

4.2 Limitations of the study

The present study while bringing insights to the impact of climate change on the AS OMZ, also possesses some unavoidable limitations. One of the main limitations arises from the simplicity of our biogeochemical model. The NPZD coupled with ROMS is based only on nitrogen and do not represent any other limiting nutrients. This can potentially lead to the biases in the representation of phytoplankton in the regions such as western AS, where nutrients other than nitrate may be a limiting factor (Naqvi et al., 2009). However, the productivity in a large part of AS is limited by nitrate (Koné et al., 2009) and hence the choice of our biogeochemical component should not affect the main conclusions of the study. We also find an underestimation of the intrusion of the Gulf water in the northern AS that introduces strong stratification there. This has already been noticed in the earlier studies with same configuration (Lachkar et al., 2018) and is attributed to inaccurate forcing fields in the Gulf region, that does not allow the correct representation of circulation and the outflow of the Gulf water into AS. However, compared to the highly biased CMIP5 models, our configuration captures the main patterns of key physical and biogeochemical patterns in the AS basin including their seasonal modulations, indicating the above factors do not necessarily limit exploiting the regional model for current attempt.

5 Conclusions

We utilize simulations from a subset of 3 selected Coupled Model Intercomparison Project Phase 5 (CMIP5) along with a set of regional downscaled experiments using ROMS configured for the Arabian Sea (AS) domain to assess the impact of climate change on the AS Oxygen Minimum Zone (OMZ), and to investigate its main drivers. The regional downscaling approach allows us to design and perform various sensitivity experiments with the same base model for the selected parent CMIP5 model. Thus, the sensitivity to the model representation of present-day conditions, as well as the contributions of local and remote forcing can be isolated. Moreover, the regional downscaling approach also allows us to correct the biases in the earth system models.

Overall, our analysis reveals a robust oxygenation in the southern AS that can potentially lead to a shrinking of the AS OMZ under the climate change. The sensitivity experiments clearly reveal the dominance of remote processes in driving this oxygenation through a combination of weaker dynamical supply and reduced biological demand. The influence of local and remote processes is highly sensitive to the layer that we consider. While the local warming dominates in the upper 100m (not shown), the role of local and remote processes shows meridional heterogeneity in the sub-surface (200-700m), with local (remote) forcing driving trends in the northern (southern) AS. Using our downscaling experiments, we also show the importance of model’s present-day biases in determining the layer of deoxygenation in the northern AS. Finally, contrasting downscaling approaches with and without correction of model biases with respect to observations revealed a strong sensitivity of future O₂ projections to model simulated current conditions. Therefore, improving the representation of present-day conditions in global models is needed to reduce uncertainties in downscaled future O₂ projections.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

All the authors contributed towards forming the objectives of current study and designing the regional downscaled experiments. PV performed the simulations and analysis, wrote the manuscript with the help of ZL and ML. All authors contributed to the article and approved the submitted version.

Funding

This research was fully funded by the research grants from Center for Prototype Climate Modeling (CPCM) and Arabian Center for Climate and Environmental ScienceS (ACCESS) at New York University Abu Dhabi (NYUAD). Both PV and ZL were supported by Tamkeen through research grant CG009 to NYUAD’s ACCESS research center.

Acknowledgments

We thank the support from the Arabian Center for Climate and Environmental ScienceS (ACCESS) at New York University Abu Dhabi (NYUAD) for providing the financial support for carrying out this research work. We thank the NYUAD for facilitating high-performance computing (HPC) resources, where our downscaling experiments and the analysis were performed. We thank Benoit Marchand, Muataz Al Barwani and the whole NYUAD HPC team for technical support. We thank the ROMS/CROCO team for giving access to the model and simulations respectively. We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in Table 1 of this paper) for producing and making available their model output. For CMIP the U.S. Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. We thank Shafer Smith (NYU, New York) for the constant encouragement and support. PV thanks Alain De Verneil and Michael Mehari (NYUAD) for fruitful discussions. PV also thanks Suresh Iyyappan (CSIR-NIO, Goa) for helping to shape the initial project proposal.

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.2023.1123739/full#supplementary-material

References