Fast local warming of sea-surface is the main factor of recent deoxygenation in the Arabian Sea

deoxygenation in the Arabian Sea Zouhair Lachkar1, Michael Mehari1, Muchamad Al Azhar1,4, Marina Lévy2, and Shafer Smith1,3 1Center for Prototype Climate Modeling, New York University Abu Dhabi, Abu Dhabi, UAE 2Sorbonne Université (CNRS/IRD/MNHN), LOCEAN-IPSL, Paris, France 3Courant Institute of Mathematical Sciences, New York University, New York, USA 4Plymouth Marine Laboratory, Plymouth, UK Correspondence: Zouhair Lachkar (zouhair.lachkar@nyu.edu)

same period. More recently, using sea glider data and historical observations, Queste et al. (2018) showed an intensification of the suboxic conditions at depth in the Gulf of Oman over the last three decades. Although there has been little work done on documenting potential deoxygenation in the Arabian marginal seas (i.e., the Red Sea and the Gulf), preliminary observations suggest ongoing deoxygenation in the Gulf with recent emergence of summertime hypoxia documented there (Al-Ansari et al., 2015;Al-Yamani and Naqvi, 2019). 5 In addition to changes in O 2 , the AS has undergone several major environmental changes over the recent decades that may intensify in the future. In particular, the AS has experienced a strong warming throughout most of the twentieth century that has accelerated since the early 1990s (Kumar et al., 2009). The warming has been particularly fast in the two Arabian marginal seas (i.e., the Red Sea and the Gulf) over the last three decades with warming rates reaching up to 0.17±0.07 • C decade −1 and 0.6±0.3 • C decade −1 in the two semi-enclosed seas, respectively (Chaidez et al., 2017;Strong et al., 2011;Al-Rashidi 10 et al., 2009). This is twice to three times faster than the global average. The AS warming has been associated with important ecological and biogeochemical changes. For instance, using satellite observations and a set of historical simulations of the northern Indian Ocean, Roxy et al. (2016) found a drop of summer productivity by up to 20% between 1950 and 2005 in the western AS that they linked to surface warming and increased stratification. Furthermore, the strong Gulf warming has been linked to recent frequent mass coral bleaching events there (Burt et al., 2019). Additionally, using a set of idealized model 15 simulations, Lachkar et al. (2019) found the warming of the Gulf to induce a reduction of the ventilation of the AS OMZ, thus causing its intensification. Important changes in the regional monsoon winds have also been reported. For instance, Goes et al. (2005) have documented an increase in summer monsoon wind intensity off Somalia between 1997 and 2005. Other studies have also suggested a climate change driven intensification of the summer monsoon winds (Wang et al., 2013) and weakening of winter monsoon winds (Parvathi et al., 2017) in the northern AS. Finally, Lachkar et al. (2018) found a tight link between 20 the strength of the monsoon winds and the depth and intensity of the AS OMZ.
While these environmental perturbations may have contributed to the observed O 2 decline over the recent decades, their potential interactions and the mechanisms through which they act to modulate O 2 remain largely poorly understood. In particular, the respective roles of recent surface warming on the one hand and summer monsoon wind intensification on the other hand are yet to be quantified. Moreover, the potential contribution of the recent fast warming of the Gulf to the declining O 2 levels 25 in the northern AS has not yet been investigated. Here we reconstruct the trends in O 2 over the period between 1982 and 2010 using a high-resolution hindcast simulation of the Indian Ocean and examine their physical and biogeochemical drivers using a set of sensitivity experiments. We show that recent deoxygenation in the northern AS has been caused essentially by surface warming, in particular in the Gulf, bringing about a reduction in the ventilation of the subsurface and intermediate layers. We also show that summer monsoon intensification enhanced oxygenation of the upper ocean south of 20 • N but has contributed to 30 deoxygenation in the northern AS. These changes are likely to have important ecological and biogeochemical consequences.

Models
The circulation model is the Regional Ocean Modeling System (ROMS)-AGRIF version (http://www.croco-ocean.org) configured for the Indian Ocean. It uses the free-surface, hydrostatic, primitive equations in a rotating environment and has a terrain-following vertical coordinate system (Shchepetkin and McWilliams, 2005). To limit diapycnal mixing errors, the diffu-5 sive component of the rotated, split, upstream-biased, 3 rd order (RSUP3) tracer advection scheme is rotated along geopotential surfaces (Marchesiello et al., 2009). The non-local K-Profile Parameterization (KPP) scheme is used to parametrize subgrid vertical mixing (Large et al., 1994). The model domain covers the full Indian Ocean from 31 • S to 31 • N and 30 • E to 120 • E with a 1/10 • horizontal resolution and 32 sigma-coordinate vertical layers with enhanced resolution near the surface. The biogeochemical model is a nitrogen-based nutrient-phytoplankton-zooplankton-detritus (NPZD) model (Gruber et al., 2006). It 10 is based on a system of ordinary differential equations representing the time-evolution of seven state variables: two nutrients (nitrate and ammonium), one phytoplankton class, one zooplankton class, two classes of detritus (small and large sizes) and a dynamic chlorophyll-to-carbon ratio. The model has a module describing the cycling of oxygen as well as a parameterization of water-column and benthic denitrification (Lachkar et al., 2016). 15 The hindcast simulation is forced with ECMWF ERA-Interim 6-hourly heat fluxes, air temperature, pressure, humidity, precipitation and winds over the period from January 1982 to December 2010. Sea surface temperature is restored to AVHRR-Pathfinder and Aqua-Modis observations and surface salinity is restored to the Simple Ocean Data Assimilation (SODA) reanalysis data using kinematic heat and freshwater flux corrections proposed by Barnier et al. (1995). Initial and lateral boundary conditions for temperature, salinity, currents and sea surface height are computed from the SODA reanalysis. The 20 initial and lateral boundary conditions for nitrate and oxygen are extracted from the World Ocean Atlas (WOA) 2013 (Garcia et al., 2013a, b). Finally, we used a monthly climatological runoff and annual river nutrient discharge from major rivers in the northern Indian Ocean derived from available observations and a global hydrological model (Ramesh et al., 1995;Dai and Trenberth, 2002;Krishna et al., 2016). We restrict our analysis to the AS region extending from 5 • S to 30 • N in latitude and from 32 • E to 78 • W in longitude (Fig S1, Supporting Information (SI)). 25 The model is spun-up using climatological forcing during 58 years and is then run for four 29-year (1982-2010) repeating forcing cycles, with the first three cycles used as a part of the spin-up period (i.e., the total duration of the spin-up phase is 145 years) and the forth cycle used for analysis (similar to the forcing protocol used in the Ocean Model Intercomparison Project, Griffies et al., 2016). The climatological run is extended for an additional 29 years to quantify the artificial trends purely driven by the model drift and contrast them to trends estimated in our hindcast run. The analysis of model drift indicates 30 a very negligible drift in salinity by the end of the climatological forcing period (Fig S2). For O 2 , the model drift decreases but remains positive (ocean oxygenation) by the end of the climatological forcing spin-up period, suggesting that the estimated 4 https://doi.org/10.5194/bg-2020-325 Preprint. Discussion started: 14 September 2020 c Author(s) 2020. CC BY 4.0 License.

Model evaluation
We use available in-situ and satellite-based observations to evaluate the performance of the model in capturing the spatial and seasonal variability in key physical and biogeochemical properties. We first evaluate the mean state by contrasting a model 5 climatology computed over the ) study period to available observation-based climatologies.
The model reproduces well the observed distributions of sea surface temperature (SST) inferred from the Advanced Very High Resolution Radiometer (AVHRR) satellite data in both winter and summer seasons (Fig 1). In particular, the model captures the intense temperature gradients prevailing in both seasons across the AS. Indeed, the model reproduces the strong temperature contrasts between the AS and its marginal seas, i.e., the Gulf and the Red Sea as well as the upwelling-driven  good agreement with that from the WOA2018 climatology (Fig 1). The AS surface circulation is similarly well represented in the model in both summer and winter seasons (Fig 2). In particular, the model simulates well the reversal of the Somali Current between summer and winter seasons as well as the prominent features of the surface circulation system including the Southern Gyre, the Great Whirl and the Southwest Monsoon Current in summer and the Northeast Monsoon Current and South Equatorial Countercurrent in winter (Schott et al., 2009).

5
The model reproduces relatively well the observed surface distribution of nitrate in both seasons, despite a tendency towards slightly overestimating surface concentrations in the oligotrophic open ocean and underestimating them in the northern AS during summer (Fig 3). The large-scale nitrate distribution at depth is also well captured by the model despite some local biases ( Fig S6, SI). The model reproduces the observed high chlorophyll-a concentrations associated with the winter and summer blooms in the northern and western AS, respectively (Fig 3). Furthermore, summer high chlorophyll concentrations 10 off the Indian west coast are also relatively well captured by the model. In contrast, the model tends to substantially overestimate chlorophyll levels off the Somali coast and in the western AS in general. This discrepancy may result from the fact that the model does not represent iron and silicic acid that may be contributing to limiting productivity in this region (Koné et al., 2009). Finally, the model captures fairly well the large-scale distribution of oxygen in the AS region (Fig 4). In particular, the location, size and intensity of the AS oxygen minimum zone is relatively well reproduced. 15 A more quantitative assessment of the model performance is shown using Taylor diagrams (Taylor, 2001) that summarize three important statistics: (i) the correlation coefficient between the model and the observations, (ii) the standard deviation in the model relative to the observations and (iii) the centered root mean square (RMS) with respect to observations (Fig 5).
To this end, in-situ observations from the World Ocean Database (2018) were binned temporally and spatially into a 0.5 • x 0.5 • seasonal climatology covering the AS region. No spatial interpolation was performed and the model and observations are 20 contrasted at the observation points only. This analysis confirms the relatively good agreement between model and observations. 8 https://doi.org/10.5194/bg-2020-325 Preprint. Discussion started: 14 September 2020 c Author(s) 2020. CC BY 4.0 License.
More specifically, it shows that surface temperature and salinity, as well as upper (400 m) ocean nitrate and O 2 distributions generally all agree well with observations with correlations above 0.9 and comparable standard deviations, across all seasons ( Fig 5). The model shows a weaker performance though in simulating surface chlorophyll with correlations ranging from 0.42 during winter months to 0.67 during the Fall season ( Fig 5).
Finally, we evaluate the model performance in reproducing observed long term changes. To this end, we binned observations 5 from the World Ocean Database (2018) seasonally on a 5 • x 5 • horizontal grid at each standard depth. We then computed seasonal climatologies for the first five years (1982)(1983)(1984)(1985)(1986) and the last five years (2006)(2007)(2008)(2009)(2010) for each bin where at least three observations were available from each year in the two periods. Using these requirements allows us to reduce the local noise in the estimated long-term changes. Yet, only temperature was found to have a coverage that is sufficiently dense to enable these anomalies to be estimated over a In summary, despite a few identified biases, the model exhibits reasonable skill in reproducing the mean seasonal climatological state for both physical and biogeochemical properties. Furthermore, it reproduces fairly well the structure and intensity of the AS oxygen minimum zone. Finally, both the magnitude and patterns of the simulated upper ocean warming are consistent 15 with observations.

Deoxygenation trends in the Arabian Sea
To investigate long-term changes in O 2 time series, data was deseasonalized by removing monthly climatologies from the original time series. Furthermore, to identify significant trends in O 2 we used the non-parametric Mann-Kendall (MK) test 20 (Mann, 1945;Kendall, 1948) that does not assume normality of data distribution and hence is less sensitive to outliers and skewed distributions. The analysis of oxygen trends between 1982 and 2010 shows a decline of O 2 concentrations in the northern and western AS in the upper 200 m (Fig. 6a). Below 200 m, the drop in O 2 is particularly important in relative terms in the northern AS but other areas in the southern and southeastern AS also show negative trends although with a weaker magnitude (Fig. 6b). A vertical transect at 65 • E indicates that in absolute terms most of the O 2 decline is concentrated north of 25 20 • N in the top 300 m (Fig. 6c). This deoxygenation results in a significant intensification of the OMZ over the three decade study period, with the volume of suboxic (O 2 < 4mmol m −3 ) water increasing by nearly 10% per decade (Fig. 7a). This causes an amplification of denitrification in the region with an increase in denitrification rates by around 14% per decade over the same period (Fig. 7c). Despite deoxygenation trends dominating in the northern and western AS, local oxygenation patches are simulated in the eastern, central and southern AS (Fig 6). This results in the net hypoxic volume (O 2 < 60 mmol m −3 ) to 30 change little over the study period (Fig 7b). Next, we focus on the northern AS (north of 20 • N) where simulated deoxygenation and its impacts on the OMZ are the most prominent.  radius (distance to the origin point) represents the modeled standard deviation relative to the standard deviation of the observations. The angle between the model point and the X-axis indicates the correlation coefficient between the model and the observations. Finally, the distance from the reference point to a given modeled field represents that field's centered root mean square (RMS) with respect to observations.
To further explore the drivers of ocean deoxygenation in the northern AS, we performed an O 2 budget analysis north of 20 • N in the 100-200 m and the 200-1000 m layers (Fig. 9). To this end, we quantified the cumulative O 2 anomalies in each of the two vertical layers and calculated the respective contributions of transport and biology to these. In order to explain changes SI). Next, we explore the changes in physical forcing causing these trends.

Impact of changes in atmospheric forcing
The analysis of long-term trends in atmospheric forcing reveals a widespread warming of the sea surface by between 0.5 and 15 1 • C in the northern AS and by up to 1.5 • C in the Gulf and in the northern part of the Red Sea, between 1982 and 2010 ( Fig   10a). In addition to surface warming, surface winds have undergone important changes with an intensification of upwelling favorable winds off Somalia and Oman, in particular during the summer monsoon season (Fig 10b and Fig S10, SI).
To quantify the contributions of surface warming and wind changes to northern AS deoxygenation, we performed three additional sensitivity experiments. In the first simulation, S hclim , all atmospheric and lateral boundary conditions were set to 20 vary interannually like in the control run except the heat fluxes that were extracted from a normal year, 1986, and repeated every year (i.e., climatological heat fluxes). In the second sensitivity run, S hclim_AG , the heat fluxes were similarly extracted from year 1986 and repeated annually, but only over the Gulf region (i.e., climatological heat fluxes over the Gulf only). Finally, in a third simulation S wclim_JJAS all atmospheric and lateral boundary conditions were set to vary interannually like in the control run except summer monsoon winds that were extracted from the year 1986 and repeated every year (i.e., climatological between 200 and 1000 m, but deoxygenation rates are nearly twice (0-200 m) to five times (200-1000 m) weaker than in the control (Fig 11 and Fig S11, SI). This indicates that the fast warming over the Gulf has likely contributed to the northern AS recent strong deoxygenation. Finally, in the absence of interannual changes in summer winds, simulated deoxygenation rates in the northern AS are weaker than in the control simulation (Fig 11). However, deoxygenation trends become stronger in most of the rest of the AS in the top 200 m and in the central and eastern AS below 200 m (Fig 11 and Fig S11, SI). This suggests 5 that summer upwelling intensification contributes to deoxygenation in the northern AS and in the western AS at depth but acts to oxygenate the upper ocean in the rest of the AS.
In summary, we conclude that deoxygenation in the northern AS is essentially caused by surface warming and that the fast warming of the Gulf plays an important role in this trend. Summer monsoon intensification contributes to deoxygenation in the northern AS as well as in the western AS at depth. However, increasing summer winds do also enhance the oxygenation of the upper ocean south of 20 • N. Next, we explore the mechanisms through which surface warming and enhanced summer monsoon winds cause these oxygen trends.

Mechanisms of Arabian sea reduced ventilation
The strong surface warming resulted in an increase in vertical stratification that is particularly important in the northern AS with 5 stratification at 200 m increasing on average by nearly 4% per decade (Fig 12a and Fig S12a, SI). This likely contributes to the reduction in the vertical ventilation in the upper ocean as revealed by our O 2 budget analysis (Fig S9, SI). Under climatological heat fluxes, upper ocean stratification does decrease in most of the AS, and in particular in the northern AS (Fig 12c and Fig   S12, SI). This confirms that most of the increase in vertical stratification in the control simulation is caused by surface warming (Fig S13a, SI     When the heat fluxes are set to be climatological over the Gulf only, the increase in vertical stratification is similar to that in the control run over most of the AS, including the region north of 20 • N (Fig S12 and Fig S12, SI). This suggests that the reduction in the ventilation of the northern AS is not entirely caused by enhanced vertical stratification and that other mechanisms are likely at play. Lachkar et al. (2019) have shown that the AS OMZ can intensify in response to strong warming of the Gulf causing a reduced outflow of the Gulf water in the northern AS. Here, we find the depth of the Gulf water has 5 shoaled locally by up to 20 m in the Gulf of Oman in the control run relative to the S hclim_AG run (Fig S14, SI).  (Fig 12b). In contrast, under climatological summer monsoon winds the themocline depth shoals almost everywhere in the AS except in the northern region where a deepening is observed (Fig 12h). This suggests that summer monsoon wind intensification cause the thermocline depth to rise in the northern AS and deepen elsewhere (Fig S15c, SI). The relatively limited role of biology in the northern AS deoxygenation is a direct consequence of the minimal change the biological productivity has experienced in the region over the study period as well as due to the negative feedback of enhanced denitrification on O 2 consumption at depth. Indeed, the biological productivity has changed little in the northern AS as enhanced stratification on the one hand and increased summer upwelling on the other hand have opposing effects on nutrient supply to the euphotic zone. For instance, productivity increases substantially in the northern AS and across the domain in the 30 absence of surface warming (Fig 13). Conversely, under climatological summer winds productivity tends to decrease in most of the northern and western AS (Fig 13). Finally, the enhanced denitrification associated with the deoxygenation trends reduces aerobic O 2 consumption and hence opposes deoxygenation in the northern AS (Fig S19, SI).

Comparison with previous works
The fact that recent subsurface deoxygenation in the northern AS has been mostly driven by reduced ventilation is consistent with previous results from Earth System models showing that tropical subsurface ocean O 2 changes projected for the twentyfirst century and those simulated during the last deglaciation are more driven by changes in ventilation than biology (Bopp et al., 2017). Our analysis has revealed that summer monsoon wind increase over the study period has mostly caused upper 5 ocean oxygenation (south of 20 • N), and to a lesser extent has contributed to deoxygenation at depth (through enhanced O 2 consumption). This is consistent with the findings of Lachkar et al. (2018) who investigated the impacts of idealized perturbations in monsoon wind intensity on the size and intensity of the AS OMZ and found that the OMZ expands and deepens in response to monsoon wind intensification. Our finding that the suboxic volume and denitrification are highly sensitive to deoxygenation is consistent with previous studies that suggest high vulnerability of suboxic zones to small changes in the .5 emission scenario . The CMIP6 multi-model ensemble average indicates an even stronger oxygen decline of more than 30 mmol m −3 by 2100, in the northern AS between 100 m and 600 m, essentially driven by an increase in the AOU . These changes can have dramatic impacts on O 2 -sensitive species and nitrogen and carbon cycling in the region. Furthermore, the impacts of these O 2 changes on marine ecosystems are likely to be exacerbated by concurrent stressors such as warming, declining productivity and ocean acidification. For instance, the fast warming of the 25 AS upper ocean is likely to cause an increase in the metabolic rates. This in turn is expected to increase organisms oxygen demand at the same time where O 2 supply to the environment is dropping (Levin, 2018;Deutsch et al., 2015). Furthermore, oxygen-induced vertical plankton migration can intensify O 2 depletion at depth (Bianchi et al., 2013). It is also established that ocean acidification interacts with deoxygenation and can aggravate the effects of low oxygen on organisms (Gobler and Baumann, 2016;Miller et al., 2016).

Caveats and limitations
Our study has a couple of caveats and limitations. Among the study's main limitations is the relatively short simulation period that precludes the attribution of the documented O 2 changes to climate change vs. natural variability. Previous observations suggest that the natural variability in O 2 , dominated by interannual and decadal oscillations, can locally be stronger than the long-term trends associated with climate warming (Whitney et al., 2007;Cummins and Ross, 2020). Strong modulation of 5 interrannual variability in hypoxic and suboxic volumes by decadal oscillations has been documented in previous studies (e.g., Deutsch et al., 2011Deutsch et al., , 2014. This complicates the detection and attribution of long-term responses to climate change. Therefore, it is not unlikely that an important fraction of the trends simulated here is associated with natural variability as it has been shown that the emergence of the climate change signal among the internal variability range is generally slow for oxygen, with only a small fraction of the ocean experiencing emergence before the end of the century (Frölicher et al., 2016). According to the 10 same study, the earliest emergence of the O 2 signal in the AS is expected to occur only by the middle of the current century.
Yet, according to other studies the emergence of the climate change signal from the noise of natural variability can occur as early as 2010 or 2020 in the ventilated thermocline of the northern AS (Long et al., 2016;Hameau et al., 2019).
Other study limitations pertain to the biogeochemical model assumptions. For instance, the lack of a representation of some major limiting nutrients such as iron, silicate and phosphate can potentially cause biases in regions where these nutrients 15 contribute to limit biological production (e.g., off Somalia). However, previous studies suggest nitrogen is the main limiting nutrient in the Indian Ocean at larger scales (Koné et al., 2009). Other major model-related limitations concern the lack of representation of important biogeochemical processes such as N 2 fixation and the crude representation of microbial respiration in the model. Yet, on the one hand recent studies suggest that N 2 fixation has a limited effect on the AS OMZ as it constitutes only a negligible proportion of new nitrogen there (Guieu et al., 2019). On the other hand, simple representations of microbial 20 respiration that miss potentially important biogeochemical feedbacks are a common problem in most existing biogeochemical models Robinson, 2019). Additional work that combine expanding the current oxygen-measurement system and improving the complexity of microbial respiration in numerical models is needed to reduce biases in model estimates of deoxygenation .

25
We reconstruct the evolution of dissolved ocean in the AS from 1982 through 2010 using a series of hindcast simulations performed with an eddy-resolving ocean biogeochemical model forced with ERA atmospheric reanalysis. We find a significant thermocline deoxygenation in the northern and western AS, with the ocean O 2 content dropping by around 2% decade −1 and 7% decade −1 in the 0-200 m and the 200-1000 m, respectively. These changes are associated with a statistically significant increase of the volume of suboxia (O 2 < 4mmol m −3 ) and denitrification by up to 30% around and 40% over the study period, 30 respectively. Using a set of sensitivity simulations we demonstrate that deoxygenation in the northern AS has been caused essentially by a widespread warming of the sea surface, in particular in the Gulf, causing a reduction in the ventilation of the subsurface and intermediate layers. Additionally, we show that a concomitant summer monsoon intensification over the study period has enhanced the ventilation and hence the oxygenation of the upper ocean south of 20 • N but has contributed to deoxygenation in the northern AS and at depth. This is because surface warming increases vertical stratification, thus reducing ventilation of the intermediate ocean, while summer monsoon wind intensification causes the thermocline depth to rise in the northern AS and deepen elsewhere, thus contributing to lowering O 2 levels in the upper 200 m in the northern AS and increasing it in the rest of the AS. Finally, wind change-driven enhanced productivity south of 20 • N contributes to deoxygenation observed 5 at depth in the western and central AS. Our findings confirm that the AS OMZ is strongly sensitive to upper-ocean warming and concurrent changes in the Indian monsoon winds. Our results also demonstrate that changes in the local climatic forcing play a key role in regional dissolved oxygen changes and hence need to be properly represented in global models to reduce uncertainties in future projections of deoxygenation.
Acknowledgements. Support for this research has come from the Center for Prototype Climate Modeling (CPCM), the New York University