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 ).

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.

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.

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

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 ).

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.

(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).

(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.

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.