The circulation model is based on a configuration of the Regional Ocean Modeling System (ROMS) AGRIF version (Shchepetkin and McWilliams, 2005) thoroughly described in Lachkar et al. (2021). The model domain covers the entire Indian Ocean at 1/10° eddy-resolving horizontal resolution and uses 32 vertical layers in terrain-following sigma-coordinates, with enhanced vertical resolution near the surface. The marine biogeochemistry and the dynamics of the lower trophic ecosystem is represented using a nitrogen-based nutrient–phytoplankton–zooplankton–detritus (NPZD) model (Gruber et al., 2006) that uses six compartments: two nutrients (nitrate and ammonium), one phytoplankton class, one zooplankton class and two classes of detritus (small and large sizes). Previous observations suggest that nitrogen is the main limiting nutrient in the Gulf, thus supporting the use of a nitrogen-based NPZD model to investigate the dynamics of Gulf hypoxia (Grasshoff, 1976; Brewer and Dyrssen, 1985; Hashimoto et al., 1998; Al-Said et al., 2018; Al-Yamani and Naqvi, 2019). Furthermore, the model represents the cycling of oxygen as well as water column and benthic denitrification (Lachkar et al., 2016). The details of the equations that govern the evolution of O₂ in the model are provided in the Supplementary Information (SI). The model simulation is forced with the ECMWF ERA-Interim 6-hourly reanalysis data. Surface salinity is restored to the Simple Ocean Data Assimilation (SODA) reanalysis data (Carton and Giese, 2008), whereas sea surface temperature (SST) is restored to Advanced Very High Resolution Radiometer (AVHRR) Pathfinder and Aqua MODIS observations using the heat flux correction method proposed by Barnier et al. (1995). Initial and lateral boundary conditions for temperature, salinity, currents and sea surface height are derived from the SODA reanalysis, whereas nitrate and oxygen initial and boundary conditions are extracted from the World Ocean Atlas (WOA) 2013 (Garcia et al., 2013a; Garcia et al., 2013b). As the used SODA reanalysis data (version 2.3) is available only up to 2010, we have restricted the simulated period to the period from January 1982 to December 2010. Finally, we restrict our analysis to the Gulf region extending from 24°N to 30°N in latitude and from 48°E to 57°E in longitude. The model is spun-up over a period of 145 years allowing a minimal drift in the upper and intermediate ocean. The details of the spin-up strategy and model drift analysis are provided in Lachkar et al. (2021).

The performance of the model in the Arabian Sea region was thoroughly assessed in Lachkar et al. (2021). In particular, the model was shown to reproduce relatively well the physical and biogeochemical mean state as well as the observed long-term changes in temperature and stratification in the Arabian Sea. The evaluation of the performance of the model in reproducing the sea surface temperature as well as the three dimensional structure of temperature, salinity and density in the Gulf is presented in detail in the SI (Figure S1–S7). Despite the scarcity of in-situ observations, this analysis reveals that the model displays a similar range of variation in temperature, salinity and density to that depicted in observations. Moreover, the simulated progression of stratification from winter to summer is comparable to that derived from observations. Additionally, the model captures relatively well the strong observed temperature gradients across the Gulf in both summer and winter. However, the model also shows some local discrepancies with respect to observations such as an underestimation of vertical stratification in both temperature and salinity with a maximum bias in summer for temperature and winter for salinity (SI). This may stem from the fact the model and observations are covering slightly different time periods or may indicate a weaker vertical stratification in the model that may lead to an overestimation of O₂ at depth and a potential underestimation of simulated bottom hypoxia.

The volume and seasonality of the Gulf outflow (through the Strait of Hormuz) is an additional key parameter we use to assess the model skill. While this metric has been used in many previous modeling studies of the Gulf, available direct observations of the Gulf outflow are very limited and come essentially from the study by Johns et al. (2003) who deployed acoustic Doppler current profiler (ADCP) moorings to measure the intensity and seasonality of the outflow across the Strait of Hormuz. In that study, the authors investigated the exchange between the Gulf and the Arabian Sea through the Strait of Hormuz over the period from December 1996 to March 1998. They found a strong and relatively steady deep outflow below 45 m of around 0.15 ± 0.03Sv (1Sv = 10⁶ m³ s^-1) that varies from 0.08 Sv in December to 0.18 Sv in March. On top of this, a weaker but strongly varying outflow of 0.06 ± 0.02 Sv is observed near the surface (above 45 m). In agreement with these observations, our model simulates an annual deep outflow of 0.13 ± 0.05 Sv that varies between 0.06 Sv in December and 0.2 Sv in May and a surface outflow of 0.05 ± 0.02 Sv (Table 1). Taken together, the model simulates a total annual outflow of 0.17± 0.06 Sv, in a good agreement with estimate of 0.21± 0.05 Sv from Johns et al. (2003). Our simulated outflow also compares well with previous estimates based on indirect measurements (Ahmad and Sultan, 1991; Pous et al., 2004) and model simulations (Chao et al., 1992; Kämpf and Sadrinasab, 2006; Yao and Johns, 2010b) (see a listing of these different estimates in Table 1). In summary, despite a few identified local biases, the model exhibits a reasonable skill at reproducing the key characteristics of the hydrography of the Gulf as known from observations. Next, we evaluate the model representation of key biological variables, chlorophyll-a and O₂.

The evaluation of the model ability to represent biological productivity is hindered by the lack of in-situ measurements of primary production in the Gulf region. Similarly, openly available in-situ measurements of Chl-a pigment, typically used as an indicator of phytoplankton biomass are also very scarce in the region. However, we contrast here our model simulated chlorophyll-a to a set of in-situ observations collected in the Qatar exclusive economic zone in the central Gulf from six research cruises conducted in April 2015, June 2015, November 2015, February 2016, April 2016, and September 2016 (Al-Naimi et al., 2017). The comparison of the surface chlorophyll concentration throughout the year between the model and the in-situ observations depicts a good agreement both in terms of the annual average and seasonality (Figure 1). The model also reproduces well the observed depth and magnitude of the deep chlorophyll maximum (DCM) in the region both in spring and summer. However, the model underestimates the magnitude of the DCM in winter and overestimates it in the autumn season (Figure 1). Interestingly, the satellite chlorophyll data from the Ocean Color Climate Change Initiative (OC-CCI) shows a strongly significant correlation with the in-situ observations in the region, despite an important systematic overestimation of in-situ measurements [Figure 1, Al-Naimi et al. (2017)]. This overestimation is likely due to the high turbidity of the Gulf waters (classified as optically complex Case II waters) and the presence of suspended sediments that alter the reliability of the chlorophyll estimates retrieved from satellites. Nevertheless, the statistically high-correlation (r = 0.795, p< 0.001) found in the region between in-situ and satellite data suggests that remotely sensed chlorophyll observations can be used to characterize the seasonality of the Gulf biological productivity (Al-Naimi et al., 2017). Therefore, to complement our model evaluation we contrast the simulated surface chlorophyll seasonal anomalies to those derived from the OC-CCI data over the entire Gulf (Figure S8, SI). In order to minimize the potential contamination of the chlorophyll signal by suspended sediments and bottom reflectance in satellite data, we restricted this comparison to the area of the Gulf that has a seafloor bathymetry deeper than 30 m. This comparison shows a good agreement between the model and the satellite chl-a data across the Gulf, both in winter and summer seasons and highlights the fact that the simulated seasonal variability in chlorophyll is consistent with that from satellite-based observations (Figure S8, SI).

Observations of oxygen in the Gulf are very scarce. Here, we compare the model simulated O₂ with observations from the extensive survey recently conducted in the Iran exclusive economic zone from summer 2018 to autumn 2019 (Saleh et al., 2021). This comparison is shown in Figure 2. Although the observations are collected outside of the period covered by the model simulation (ending in December 2010), this comparison is useful to illustrate the model performance in reproducing the seasonal variability in O₂ both near the surface and close to the bottom. This comparison reveals that the model is able to capture the seasonal progression of low O₂ and bottom hypoxia as inferred from the Saleh et al. (2021) observations. Indeed, throughout the sampled months oxygen is lowest in the late summer and early autumn (September) below 50 m and highest near the surface in May both in the model and the observations (Figure 2). However, the model simulated O₂ concentrations are slightly above those measured in the summer and autumn, particularly below 50 m. Although not statistically significant (as the error bars from observations and model do overlap), this difference may also reflect the fact the model and observations are covering different time periods or may indicate an overestimation of O₂ at depth that may cause an underestimation of the extent and magnitude of Gulf bottom hypoxia in the model.

Finally, the evaluation of the model’s ability to reproduce long-term changes in the Gulf is limited by the sparse data coverage in the region. Nevertheless, we use available observational data to assess the model simulated trends. More specifically, we contrast here simulated trends in SST to trends based on different observational products. Additionally, we also compare simulated long-term trends in temperature, salinity and density to trends derived from the ECMWF Ocean Reanalysis System 5 (ORAS5) reanalysis at both surface and at 50 m depth. This comparison reveals that the simulated surface warming is generally in good agreement with that from the different observational SST products, despite differences in the magnitude of warming, with ERA5 (HadlSST) displaying the strongest (weakest) rates of warming (Figure S9, SI). The model comparison with ORAS5 reveals that the simulated trends in temperature at both the surface and bottom compare relatively well with those derived from the reanalysis (Figure S10, SI). However, our model simulates a slight increase in salinity, especially in the southern Gulf at the surface, that contrasts with the a slight but statistically insignificant decrease in ORAS5 salinity. This is because our model is forced by the SODA reanalysis that similarly shows a slight increase in surface salinity over the study period. This discrepancy is likely to have limited implications for the simulated trends in circulation as the water density generally decreases both near the surface and at depth in both our model and ORAS5 reanalysis (Figure S12, SI).

In addition to strong warming, the analysis of trends in atmospheric conditions over the Gulf region reveals an important weakening of the northwesterly Shamal winds that is particularly important during winter (Figure S13, SI). These trends are relatively robust and consistent across multiple data products. Indeed, trends in surface winds from three reanalysis products: ERA-Interim (here used to force the model), the Japanese 55-year Reanalysis Project (JRA-55), and the National Centers for Environmental Prediction reanalysis II (NCEP-2) all show a weakening of wind speed and a reduction in the northwesterly Shamal winds that is particularly important in winter (Figure S13, SI). Finally, a previous analysis of simulated oxygen trends in the Arabian Sea conducted over the same period revealed a strong intensification of the oxygen minimum zone (OMZ), particularly in the adjacent Sea of Oman (Lachkar et al., 2021). Therefore, to disentangle and quantify the potential individual contributions of: (i) local warming, (ii) wind changes and (iii) remote deoxygenation to the simulated O₂ trends in the Gulf, we analyze three additional sensitivity simulations (Table S1, SI). In the first run, S_wclim, all atmospheric conditions vary interannually like in the control run except winds that are extracted from a normal year (that is neutral with respect to the leading interannual climate oscillations El Niño Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD)), 1986, and repeated every year (i.e., no wind changes). This approach allows filtering out interannual variability while conserving the high-frequency variability of the original forcing (e.g., Large and Yeager, 2004). In the second sensitivity experiment, S_{AG_hclim}, the heat fluxes were similarly extracted from the year 1986 and repeated annually but only over the Gulf region (i.e., no local warming over the Gulf). Finally, in a third simulation S_{AS_hclim}, the heat fluxes are interannual over the Gulf and climatological over the rest of the domain (i.e., warming over the Gulf only). In this latest run, no significant OMZ intensification is simulated in the upper layers of the Arabian Sea (i.e., no remote deoxygenation affecting the Gulf inflow source waters).

O₂ in the control run remains close to saturation in the upper 50 m and shows a relatively limited seasonal variation, with a peak in March and a minimum in October (Figure 3A). Below 50 m, O₂ shows a much larger deviation from its saturation levels and exhibits a stronger seasonal variation with a peak in February and a minimum in October (Figure 3B). O₂ depletion at depth leads to the emergence of near-bottom (here defined as the model deepest layer) hypoxia (O₂< 60 mmol m^-3) that develops over a significant portion of the deeper part of the Gulf (bathymetry deeper than 50 m) with a maximum covered area fraction reaching up to 30-50% in the autumn (October-November) (Figure 3D). This seasonal progression of bottom hypoxia is consistent with recent observations from Saleh et al. (2021) that also indicate that near-bottom hypoxia emerges in the deepest part of the Gulf in summer and reaches its maximum spatial extent in the autumn. Our finding is also consistent with the observations of Al-Ansari et al. (2015) who first reported hypoxia in the deeper (>50 m) layers of the exclusive economic zone of Qatar in the central Gulf in the late summer and early autumn (September-October). The drop in near-bottom O₂ concentrations also leads to the emergence of suboxia (O₂< 5 mmol m^-3) during the autumn over 2-5% of the deep Gulf area (Figure 3D).

Al-Ansari et al. (2015) have speculated that bottom O₂ concentration is highest in the winter because of the weak stratification that favors the mixing of deep layers with the O₂ saturated surface waters. They also speculated that bottom hypoxia develops in summer because of the stronger stratification that prevents mixing as well as the strong fertilization of the Gulf by the intrusion of the surface nutrient-rich Indian waters in summer, that enhances productivity and boosts remineralization and O₂consumption. Here we perform an O₂ budget analysis in the deeper (>50 m) part of the Gulf to elucidate the drivers of O₂ depletion and the factors controlling the seasonal progression of hypoxia (Figure 3C). The results of this analysis indicate that near-bottom O₂ is, indeed, maximum in winter because of a weaker stratification favoring a more vigorous vertical mixing with surface waters. However, our analysis also shows that the drop in deep Gulf water oxygen levels is not driven by a peak in O₂ consumption in summer but is rather driven by the cumulative year-long biological drawdown of O₂, in combination with the weak oxygen replenishment from spring to autumn due to enhanced stratification during this period of the year (Figure 3C). Our analysis also reveals that the peak fertilization of the Gulf occurs in winter (Figure 10 and Figure S14, SI). This is likely due to the higher nitrogen content of the Indian Ocean surface waters in winter that results from convective mixing in the northern Arabian Sea (Naqvi et al., 2002; Al-Yamani and Naqvi, 2019).

The analysis of O₂ trends over the study period reveals a steady decline in the O₂ content of the Gulf as a whole as well as in its deeper (> 50 m) layer at statistically significant rates of 1.04 mmol m^-3 decade^-1 and 4.07 mmol m^-3 decade^-1, respectively (Figure 4A). The decline is particularly important in summer and autumn where O₂ concentrations decrease at depth by more than 20 mmol m^-3 between the early 1980s and the late 2000s (Figure S15, SI). The rate of O₂ decline is however uneven, with a stronger decline between 1990s and 2000s than between the 1980s and 1990s (Figure 3B), reflecting a strong imprint of the interannual variability on the long-term deoxygenation signal (Figure 4). The O₂ depletion results in an important increase in the near-bottom hypoxia- and suboxia-prone areas between the 1980s and the 2000s at statistically significant rates of 3818 and 1013 km² decade^-1, respectively (Figure 4B). The expansion of hypoxia comes with a strong seasonality with a particularly fast intensification in late summer and in the autumn seasons (Figure 3D). Indeed, not only do hypoxia (and suboxia) reach their maximum extent in the late summer and autumn but the maximum decline in O₂ in absolute terms is also recorded during this period of the year (Figure 3D and Figure S15, SI). However, in relative terms the increase in hypoxic area is more important in spring and early summer (Figure 3D and Figure S16, SI).

These changes are accompanied by a local lengthening of the hypoxic season by 1 to 3 months mostly associated with an earlier start of the hypoxic season by 1 to 2 months between the early 1980s and the late 2000s (Figure 5). Additionally, an increase in the frequency and severity of the hypoxic episodes can also be observed over the study period (Figure 5). For instance, in the late 2000s O₂ levels dropped to the suboxic or near suboxic ranges (O₂< 20 mmol m^-3) in areas that only recorded moderate hypoxia (40 mmol m^-3< O₂< 60 mmol m^-3) in the early 1980s (Figure 5A and Figure S17, SI). This has led to a near doubling of denitrification rate in the Gulf over the study period (Figure S18, SI). Additionally, in the late 2000s hypoxia became regular (with at least 1 event per year) in areas where its occurrence was uncommon (e.g., less than 1 event in 5 years) in the early 1980s (Figure 5). Next, we explore the drivers of oxygen changes in the Gulf with a focus on changes in near-bottom O₂ and hypoxia.

At the sea surface, dissolved oxygen is generally at or near-equilibrium with the atmosphere. This equilibrium concentration, O 2 sat , known as oxygen saturation level is a function of temperature, salinity and atmospheric pressure (Garcia and Gordon, 1992). As water leaves the surface and is subducted into deeper layers, its O₂ content gets depleted by biological activity. The O₂ deviation for a water parcel from its saturation level, known as apparent oxygen utilization ( AOU = O 2 sat – O 2 ), is a measure of oxygen consumption via biological activity since the parcel was last at equilibrium in contact with the atmosphere. Therefore AOU depends both on biological activity as well as on ventilation (pathways and timescales). Here, to separate the contributions of changes in ventilation and biology (O₂ utilization effect) from those driven by changes in O₂ saturation (thermal effect), essentially induced by temperature changes, we decompose the oxygen anomaly ΔO₂ as:

where Δ O 2 sat and ΔAOU represent the anomalies of O 2 sat and AOU, respectively. The analysis of the three terms of this equation reveals that the reduction of O 2 sat (thermal effect) explains the majority of the simulated O₂ decline throughout the water column in the shallow banks of the southern Gulf (Figure 6). In the western and northwestern Gulf, warming-induced decline in O 2 sat explains around half of the O₂ decline. In contrast, in the deeper part of the Gulf along the Iranian coast, changes in O 2 sat explain only 30 to 40% of the total O₂ decline, driven mostly by changes in O₂ biological utilization (AOU) (Figure 6). When considering the deepest layer of the Gulf (deepest model sigma layer), the dominance of AOU over O 2 sat is even more pronounced, with AOU changes explaining locally between 70% and 90% of O₂ decline in the deep part of the Gulf (h > 50 m) (Figure 7). We conclude therefore that the near-bottom hypoxia expansion reported in our study is mostly driven by changes in AOU, and hence changes in either biological consumption or O₂ transport, with only a small contribution emanating from decreasing O₂ solubility.

To disentangle the roles of biology (remineralization) and transport (ventilation) in driving near-bottom O₂ depletion, we performed an O₂ budget in the model deepest layer (Figure 7). This analysis reveals that near-bottom O₂ decline is driven by biological consumption, while O₂ transport acts to partially counteract this tendency (Figure 8). This compensation reflects the effect of vertical mixing with O₂-rich surface waters, dampening O₂ depletion. The O₂ budget also reveals that the enhanced O₂ consumption near the bottom is mostly driven by an increase in remineralization in the sediment (Figure 8). This reflects the importance of benthic respiration (and sedimentary processes in the Gulf biogeochemistry in general) as can be expected in such a very shallow marginal sea (Al-Yamani and Naqvi, 2019). This is also consistent with the recent findings of De Verneil et al. (2021) that suggest that a considerable fraction of the total O₂ consumption in shallow portions of the Gulf may be associated with benthic respiration.

The enhanced respiration in the Gulf is a direct consequence of productivity increase (locally reaching up to 20% of the mean production levels) over the study period (Figure 9). This increase in productivity can result from a weakening of the nutrient limitation or an increase in the temperature-dependent and light-limited phytoplankton growth rate. Examining these potential contributions, we find that most of the productivity enhancement reported here is associated with a stronger nutrient supply (fertilization) while the increase in phytoplankton growth rate (in response to water temperature increase) is playing a significant yet smaller role (Figure 9). The Gulf fertilization is indicative of an enhanced inflow of nutrients through the Strait of Hormuz as no other atmospheric or terrestrial sources of nutrients is considered in the study. To elucidate the mechanisms leading to such changes, we examine the evolution of total inorganic nitrogen (TIN) transport across the Strait of Hormuz as well as the long-term change in the Gulf water inflow (Figure 10). This analysis confirms that the transport of TIN into the Gulf substantially increased between the 1980s and the 2000s in the winter season during which nitrogen transport is largest. This is associated with enhanced surface inflow into the Gulf in winter (Figure 10). We hypothesize that this increase in the Gulf inflow is associated with the concurrent weakening of the winter northwesterly Shamal winds (Figure S13, SI) as a strong link between the Shamal winds and the strength of the cyclonic gyre in the southern Gulf (and hence the transport across the Strait of Hormuz) has been reported in previous studies (e.g., Thoppil and Hogan, 2010a; Pous et al., 2013). To test this hypothesis, we explore the sensitivity of the Gulf inflow of water and nutrients to changes in atmospheric conditions including changes in surface winds.

In order to unravel the role of (i) wind changes, (ii) local warming and (iii) remote warming and O₂ changes in the Arabian Sea in driving the Gulf O₂ decline, we contrast the control run with three additional sensitivity simulations. More specifically, we identify the effect of wind changes on oxygen by contrasting the changes in the no-wind-change S_wclim run to those in the control run. Similarly, contrasting the O₂ changes in the no-Arabian Sea-warming simulation S_{AS_hclim} to those in the control run allows us to measure the relative importance of the warming of the Arabian Sea to Gulf deoxygenation. As O₂ changes in the northern Arabian Sea are non-significant in the S_{AS_hclim} run [see Lachkar et al. (2021)], this comparison also reveals the importance (or the lack thereof) of deoxygenation trends in the Arabian Sea to the O₂ decline in the Gulf. Finally, by comparing oxygen changes in the S_{AG_hclim} (no Gulf warming) and the control run, we are able to measure the effect of local Gulf warming on oxygen changes and hypoxia expansion. This analysis confirms that changes in winds (weakening of Shamal winds) are responsible for most of the O₂ decline at depth in the Gulf (Figures 11, 12), and therefore are the primary cause of near-bottom hypoxia expansion (Figure 13). Indeed, under climatological winds (i.e., S_wclim) near-bottom O₂ depletion is much weaker than in the control in both summer and autumn seasons (Figure 11). Additionally, no significant trends can be identified in either O₂ content or hypoxic area in the absence of wind changes (Figures 12,  13). Contrasting O₂ changes in the no-Arabian Sea-warming simulation S_{AS_hclim} to those in the control run reveals very limited differences, except around the strait of Hormuz region and in the eastern Gulf where there is a weaker decline in O₂ relative to the control (Figures 11,  12). The expansion of the hypoxic area is also only marginally weaker in the S_{AS_hclim} run relative to the control (Figure 13). This indicates that while the warming of the Arabian Sea (and its deoxygenation) may slightly enhance O₂ depletion in the eastern Gulf, it has a very limited impact on deoxygenation in the rest of the Gulf. The comparison of the S_{AG_hclim} (no Gulf warming) and the control run reveals important differences near the sea surface with weaker O₂ decrease in the upper 20 m in the absence of Gulf warming (Figure 11). However, O₂ decline in the deep layers is comparable in the two simulations (Figures 11, 12). Furthermore, near-bottom hypoxia expansion is relatively similar in both simulations (Figure 13). We therefore conclude that while local warming of the Gulf leads to O₂ loss in the near-surface due to reduced O₂ solubility and outgassing, it has a relatively limited impact on deep O₂ levels and contributes little to simulated bottom hypoxia expansion. Finally, contrasting N O 3 − changes simulated in the control run to those obtained under different atmospheric forcing scenarios further confirms that the Gulf fertilization is linked to surface wind changes (Figure S19). Indeed, in the absence of wind changes (i.e., S_wclim), no increase in N O 3 − is observed in the Gulf (Figure S19). Additionally, in contrast to the control run and the sensitivity runs where surface winds do change (i.e., S_{AS_hclim} and SA_{G_hclim}) that show an increase in the surface water inflow through the Strait of Hormuz over the study period, no similar change is produced in the absence of wind changes (Figure S19).

In summary, we simulate an important deoxygenation in the Gulf that causes the intensification and expansion of the seasonal near-bottom hypoxia in the deeper layers. Our analysis reveals that O₂ decline is mostly driven by local warming near the sea surface and by the weakening of Shamal winds that enhances the Gulf fertilization in the winter. This results in increased production in the surface and enhanced respiration near and in the sediment, causing near-bottom oxygen depletion (Figure 14). Deoxygenation in the Arabian Sea appears to have limited impact on Gulf oxygen content, beyond amplifying O₂ depletion in the eastern Gulf.

Observational evidence suggests an expansion and intensification of near-bottom hypoxia in the Gulf (El Samra and El Gindy, 1990; Al-Ansari et al., 2015; Al-Yamani and Naqvi, 2019; Saleh et al., 2021). Al-Ansari et al. (2015) have estimated the area covered by hypoxic water in the central Gulf in summer to exceed 7,000 km². However, this is likely an underestimation as this study was restricted to the Qatar exclusive economic zone. Saleh et al. (2021) have observed hypoxic waters in the western Gulf (November 2018) at depths exceeding 50 m and in the eastern Gulf at depths exceeding 75 m (December 2019). Assuming that hypoxia covered the entire seafloor of the Gulf at depths below 50 m in the western part and 75 m in the eastern part, these authors estimated the total area of the Gulf prone to bottom hypoxia in mid-autumn to potentially exceed 50,000 km². Here, we report an increase of the hypoxic area from less than 20,000 km² in the early 1980s to above 30,000 km² in the 2000s. (Figure 4B). The discrepancy between the two hypoxia-prone area estimates may reflect an excessive extrapolation by Saleh et al. (2021) or stem from an underestimated bottom hypoxia in our model, potentially caused by the simulated weaker vertical stratification relative to observations.

Both Al-Ansari et al. (2015) and Saleh et al. (2021) reported high phosphate and nitrate concentrations in the bottom hypoxic waters, suggesting that bottom hypoxia is tightly related to enhanced remineralization in the benthos and near the seafloor. De Verneil et al. (2021) also reported that a substantial fraction of O₂ consumption at a reef in the southern Gulf is associated with benthic respiration. These results are in agreement with our finding that O₂ depletion near the bottom is mostly associated with the increase of remineralization in the sediment.

Al-Ansari et al. (2015) have speculated that summer hypoxia is driven by the fertilization of the Gulf in summer by the inflow of Indian Ocean surface water driven by summer monsoon winds. Our study shows the summer and autumn hypoxia in the Gulf is caused by the strong stratification prevailing during this period of the year, preventing O₂ replenishment of the bottom waters combined with the persistent O₂ consumption driven by year-long remineralization near the bottom and in the benthos. Additionally, our analysis also suggests that the Gulf fertilization is stronger in winter than summer because of the higher nutrient content of the Indian Ocean surface waters in winter. This is consistent with observations suggesting winter convective mixing in the northern Arabian Sea and the Sea of Oman to cause nutrient entrainment into the surface layer (Madhupratap et al., 1996; Naqvi et al., 2002). In contrast, summer Indian Ocean surface waters are more nutrient-depleted because of the stronger stratification and because upwelled waters along the southern coast of Oman do not enter the Sea of Oman (Naqvi et al., 2006).

Al-Said et al. (2018) reported elevated levels of total organic carbon (TOC) in the Kuwaiti waters between 2014 and 2016 linked to anthropogenic pollution. These authors have suggested that the recent expansion of hypoxia in the Gulf might be driven by increased O₂ biological consumption associated with eutrophication. Our analysis also indicates that the intensification of the O₂ depletion in the Gulf is primarily driven by enhanced respiration at depth. However, our results show that enhanced respiration in the Gulf is not necessarily driven by terrestrial pollution but can also be caused by larger scale changes in climatic factors.

Finally, using a combination of in-situ observations and modeling, Paparella et al. (2019) and Burt et al. (2019) have demonstrated that summer bottom temperature in the Gulf is tightly controlled by the intensity of evaporative cooling driven by Shamal winds. For instance, the two studies linked the dramatic warming of the sea bottom recorded in summer 2017 (that led to a severe coral bleaching event in the southern Gulf) to the absence of strong Shamal wind events in that particular summer. Our study sheds light on the role of weaker winter Shamal winds in the enhanced fertilization of the Gulf and hence the intensity of subsequent biological respiration and O₂ depletion. A weakening of the summer Shamal winds as suggested by Figure S13 (SI) can contribute to further amplify O₂ depletion at depth, as the lack of wind-driven evaporative cooling is likely to increase surface warming and stratification in the deeper part of the Gulf (causing a weaker O₂ replenishment at depth) as well as lead to an increase in bottom temperature in the shallower Gulf (causing a decreased O₂ saturation).

A recent study by De Verneil et al. (2021) reported the occurrence of bottom hypoxia at very shallow depths (∼ 4 m) off the coast of Abu Dhabi in the summer and early autumn of 2019. However, this consisted in short-lived (median 3 hours) night-time hypoxic episodes caused by overnight respiration consuming up most of the available dissolved O₂. Our study shows that in the shallow portion of the Gulf, long-term changes in O₂ are equally sensitive to declining O₂ solubility (and increasing water temperatures) as to changes in O₂ biological utilization. Therefore, it can be expected that further warming of the Gulf in the future may lead to an increase in the frequency and possibly the duration of these short-lived but recurrent hypoxic events, which is likely to have strong implications for marine organisms and sensitive ecosystems such as coral reefs, as well as for economically important sectors such as fisheries (De Verneil et al., 2021; Maltby et al., 2022).

The hypersaline dense Gulf waters sink to intermediate depths (200-300 m) when they enter the Arabian Sea and hence are essential to the ventilation of its OMZ (McCreary et al., 2013; Lachkar et al., 2019). It has been recently shown that the warming of the Gulf waters limits this ventilation, as the warmer and more buoyant Gulf waters spread at shallower depths in the Arabian Sea, thus leading to an intensification of the upper part of the OMZ (Lachkar et al., 2019). Al-Said et al. (2018) and Al-Yamani and Naqvi (2019) have speculated that an increase in the TOC content of the Gulf waters, driven by increasing levels of anthropogenic pollution or eutrophication in the Gulf, can ultimately lead to an intensification of the Arabian Sea OMZ. This is because the semi-labile fraction of the TOC can in principle survive during the transit of water in the Gulf and become available for decomposition once the Gulf waters exit the Hormuz Strait and are exported into the OMZ. Here we argue that declining O₂ concentrations in the Gulf, driven by changes in the Gulf climate as shown here, may further amplify the intensification of the Arabian Sea OMZ. This is a potential additional mechanism (in addition to (i) the impact of Gulf warming on the Arabian Sea ventilation and (ii) the impact of pollution-driven increase in TOC levels) through which local changes in the Gulf physical and biogeochemical environment can affect the Arabian Sea OMZ.

O₂ decline is only one among many stressors that the Gulf ecosystems have to deal with. Rapid warming and additional changes in pH (acidification), salinity (desalination) and nutrient load (eutrophication) are a few other examples of stressors that may further push these ecosystems closer to critical environmental thresholds (Le Quesne et al., 2021; Lincoln et al., 2021; Maltby et al., 2022). Moreover, these multiple stressors may interact in complex ways and lead to synergistic or antagonistic effects that cause a larger or a smaller impact on the organisms and ecosystems than when experienced in isolation. Understanding these interactions is key to narrowing down uncertainties around the impact of future climate change on the Gulf ecosystems and to devising effective mitigation strategies for conservation of sensitive species threatened by oxygen decline.

Finally, our study reveals that hypoxia in the Gulf is not only subject to strong seasonality but is also associated with a strong interannual variability. Therefore, improving our understanding of the dynamics of dissolved O₂ in the Gulf and documenting its variability requires devising suitable observational strategies that help ensure a sufficient sampling in order to detect long-term changes. This can be achieved through deploying fixed oceanographic moorings that establish long time-series and ensure a continuous physical and biogeochemical monitoring of the Gulf region following successful examples from other subtropical systems such as the Bermuda Atlantic Time-series Study (BATS) in the western subtropical Atlantic or the European Station for Time series in the Ocean Canary Islands (ESTOC) in the subtropical Eastern Atlantic. Alternatively, devising an observational program that consists of collecting key observations along repeated transects at regular intervals of time and ideally during the hypoxic season (summer-autumn), may also provide a wealth of information on the processes responsible for O₂ depletion in the Gulf, but also on the characteristics of hypoxia (e.g., intensity, frequency, duration, etc.) and its biogeochemical and ecological implications.

Our study has a couple of caveats which include the limitations associated with the modeling setup. For instance, while the model horizontal resolution (1/10°deg) is 4 to 5 times smaller than the local Rossby deformation radius (Chelton et al., 1998), allowing to explicitly resolve most mesoscale eddies, it remains relatively coarse to represent the fine details of the geometry of the coastline and the seafloor bathymetry. Although previous modeling studies at comparable resolutions (e.g., Kämpf and Sadrinasab, 2006; Pous et al., 2015) have been shown to reproduce the essential characteristics of the circulation and hydrography of the Gulf, enhanced resolution may be needed to improve the realism of the local currents and the influence of mesoscale eddies on the mean circulation and mixing as well as to refine the representation of stratification (Pous et al., 2015; Lévy et al., 2022). An additional model related caveat pertains to the absence of external sources of nutrients either due to coastal pollution and eutrophication (Devlin et al., 2015) or via atmospheric dust deposition (Al-Yamani and Naqvi, 2019). The crude representation of microbial respiration in the model is another major model-related limitation. For instance, remineralization rates in the model are set to be constant, while laboratory studies suggest that the metabolism processes that are responsible for respiration are temperature-dependent (López-Urrutia et al., 2006). Although simplistic representations of microbial respiration are a common problem in most existing biogeochemical models (Laufkötter et al., 2017; Oschlies et al., 2018), additional work to enhance the realism of this key process in the model is needed to help reduce uncertainty around model simulated O₂ trends in the region.

Second, our analysis reveals that the evolution of O₂ in the Gulf is marked by a strong interannual and decadal variability. Therefore, the trends reported in this study may not necessarily be representative of the long term climate-change driven perturbation. Previous observations suggest that the natural variability in O₂, dominated by interannual and decadal oscillations, can locally be stronger than the long-term trends associated with climate warming (e.g., Whitney et al., 2007). This is particularly true as productivity and O₂ consumption in the Gulf are very sensitive to changes in water exchange through the Strait of Hormuz, which in turn is strongly sensitive to interannual variability in atmospheric fluxes (Pous et al., 2015).

The prescribed main changes in the atmospheric conditions, namely the Gulf warming and changes in Shamal winds are relatively robust features across different data-based products (Figure S9 and S13, S1). However, the simulated changes in productivity and nutrient supply to the Gulf remain to be confirmed by in-situ observations. In particular, repeated measurements of key physical parameters such as the intensity and structure of the water exchange through the Strait of Hormuz and important biological parameters such as chlorophyll and nutrient concentrations are essential to documenting ongoing changes and confirming the mechanisms described in our modeling study.