Effects of upwelling duration and phytoplankton growth regime on dissolved-oxygen levels in an idealized Iberian Peninsula upwelling system

. We apply a coupled modelling system com-posed of a state-of-the-art hydrodynamical model and a low-complexity biogeochemical model to an idealized Iberian Peninsula upwelling system to identify the main drivers of dissolved-oxygen variability and to study its response to changes in the duration of the upwelling season and in the phytoplankton growth regime. We ﬁnd that the export of oxygenated waters by upwelling front turbulence is a major sink for nearshore dissolved oxygen. In our simulations of summer upwelling, when the phytoplankton population is generally dominated by diatoms whose growth is boosted by nutrient input, net primary production and air–sea exchange com-pensate dissolved-oxygen depletion by offshore export over the shelf. A shorter upwelling duration causes a relaxation of upwelling winds and a decrease in offshore export, resulting in a slight increase of net dissolved-oxygen enrichment in the coastal region as compared to longer upwelling


Introduction
Declining dissolved-oxygen (DO) levels in the world's oceans, including in the coastal realm, impact all marine life (Laffoley and Baxter, 2019), from microbes to higher trophic levels (Breitburg et al., 2018), with consequences ranging from ecological adaptations and shifts (Gilly et al., 2013) and changes of biogeochemical activity (Wright et al., 2012) to mass mortality events (Diaz and Rosenberg, 2008) and biodiversity restructuring (Vaquer-Sunyer and Duarte, 2008).
Coastal waters are generally eutrophic and characterized by substantial planktonic productivity at the surface which favours oxygen consumption through the remineralization of sinking organic matter, leading to low levels of DO in subsurface and near-bottom waters. Highly productive surface waters are typical of coastal upwelling regions and sustain socioeconomically important ecosystems. Upwelling systems are thus especially sensitive to episodic and longterm changes in DO levels with deleterious effects on marine life and human activities such as fisheries (Grantham et al., 2004;Hales et al., 2006;Paulmier and Ruiz-Pino, 2009;Mc-Clatchie et al., 2010;Roegner et al., 2011).
The western Iberian Peninsula Upwelling System (IPUS) is the northern branch of the Canary Upwelling System where the intra-annual variability of alongshore winds produces a seasonal upwelling-downwelling cycle (Wooster et al., 1976). In the summer and early fall, the Azores High migrates northward, causing a poleward alongshore wind that forces offshore Ekman transport of surface waters and the upwelling of subsurface waters, while downwelling prevails the Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
rest of the year. Long-term studies of the seasonal upwelling in the IPUS have pointed to a weakening of upwelling winds over multi-decadal timescales (Sousa et al., 2017), but simulations of future climate indicate an enhancement of the upwelling due to the poleward migration of the Azores High and a lengthening of the upwelling season (Miranda et al., 2013;Rykaczewski et al., 2015;Sousa et al., 2017).
During the upwelling season, hydrographic and biogeochemical variability is primarily determined by the wind forcing that controls the inflow of offshore subsurface water masses onto the shelf (Álvarez-Salgado et al., 1993). DO levels in these subsurface waters are thought to diminish during the upwelling season due to the strong remineralization of sinking organic matter as well as to mixing with poorly ventilated waters of subtropical origins (Castro et al., 2006).
Overall, DO appears to partly control the remineralization of dissolved organic carbon in the IPUS and more generally in key water masses of the North Atlantic (Álvarez-Salgado et al., 2013).
DO, in addition to being transported by the Ekman upwelling circulation, is also affected by the turbulent component of the circulation, characterized by (sub-)mesoscale fronts, filaments and eddies (Bettencourt et al., 2015;Capet et al., 2008;Chaigneau et al., 2009;Marchesiello et al., 2003;Montes et al., 2014;Vergara et al., 2016). These structures have a strong influence on the cross-shore transport and on the vertical redistribution of biogeochemical tracers (Bettencourt et al., 2017;Combes et al., 2013;Gruber et al., 2011;Hernández-Carrasco et al., 2014;Nagai et al., 2015;Renault et al., 2016;Rossi et al., 2013). In the Iberian basin, the Eastern North Atlantic Central Water (ENACW) mass of subtropical origin, whose archetypal concentrations of DO are initially low, would be ventilated by eddy-induced mixing with more oxygenated ENACW of subpolar origin (Pèrez et al., 2001).
From the biological perspective, the planktonic community structure can also affect DO in the water column, as this results from the balance between oxygen production by photo-autotrophs and respiration by both autotrophs and heterotrophs. Changes in the community composition should lead to changes in the DO budget of the continental shelf. In the Iberian upwelling, planktonic blooms tend to be dominated by large cells (micro-phytoplankton), such as diatoms (Moncoiffé et al., 2000;Rossi et al., 2013), which have lower respiration-to-photosynthesis ratios than dinoflagellates or cyanobacteria (López-Sandoval et al., 2014). Moreover, diatoms growth is thought to be positively correlated with the upwelling-driven nitrate inputs (Sarthou et al., 2005); in contrast, the growth of other planktonic groups can be insensitive to and even limited by newly upwelled nitrates (Mahaffey, 2005), especially when colimitation prevails or in the absence of necessary micro-nutrients. Thus, changes in the relative dominance of functional groups in the community can lead to different autotrophy-heterotrophy regimes and total DO levels.
In this paper, we perform a process-oriented study coupling an idealized IPUS configuration of a hydrodynamic model with a low-complexity biogeochemical model centred on oxygen. We focus on the mechanisms of DO change to assess the relative role and importance of physical (upwelling duration) and biogeochemical (functional planktonic community) processes in setting DO levels and distribution. Previous studies have focused on the seasonal and climatological modelling of the planktonic ecosystems and of dissolved oxygen (Marta-Almeida et al., 2012;Reboreda et al., 2014Reboreda et al., , 2015. Here, the spatiotemporal scales studied are shorter, since we are interested in the variations and forcing mechanisms of DO levels in response to successive and short-lived upwelling pulses, commonly observed in the IPUS during the upwelling season. The coupled physical-biogeochemical model is presented in Sect. 2; our results are shown and analysed in Sect. 3. We then discuss our findings and conclude in Sect. 4.

The coupled model
The idealized IPUS is contained in a meridionally oriented stretch of continental shelf with constant bathymetry and limited zonally by the coast and the open ocean (Fig. 1). The domain has 300 km of length and width, with a shelf width of 50 km, which is characteristic of the Iberian shelf at 41 • N. The resolution of the spatial grid is 500 m, allowing for the simulation of summer upwelling seasons up to 90 d duration, with the associated upwelling front dynamics -mesoscale and submesoscale vortices and filaments -within acceptable computing time limits.
The physical-biogeochemical model couples the 3D Regional Ocean Modelling System (ROMS; Shchepetkin, 2015;Shchepetkin and McWilliams, 2005), in its Coastal and Regional Ocean Community (CROCO) version (Debreu et al., 2012), to the low-complexity oxygenphytoplankton-zooplankton biogeochemical model of Petrovskii et al. (2017) and Sekerci and Petrovskii (2015), hereafter SP2015. Oxygen is used as the main model currency, since it is central for the functioning of marine ecosystems (Breitburg et al., 1997(Breitburg et al., , 2018. Oxygen is also at the heart of the photosynthesis and carbon fixation metabolic pathways for prokaryotic and eukaryotic algae, being a major electron acceptor (Mackey et al., 2008;Zehr and Kudela, 2009) and is needed for many biogeochemical reactions, including those involved in remineralization, which constantly occurs in the oceanic water column. The low-complexity O 2 PZ model was directly implemented into the CROCO model following one of the already-embedded biogeochemical modules; it allows for maintaining the full capabilities and versatility of this community code while gaining in computing efficiency (allowing for extensive sensitivity studies and the consideration of submesoscale dynamics in oxygen cycling) as compared to the more complex biogeochemical modules already available within CROCO.
ROMS-CROCO is a hydrostatic primitive equation model, formulated in finite differences with time splitting in the barotropic (fast) mode and the baroclinic (slow) mode. It uses a terrain following vertical discretization in σ coordinates with 40 vertical levels. The prognostic variables are the horizontal barotropic and baroclinic momentum components, temperature and salinity. Free-surface displacement and vertical momentum are diagnostic variables. The horizontal advection of momentum and tracers is calculated from a third-order upstream-biased scheme. Horizontal mixing is Laplacian, along iso-σ surfaces, with diffusivity coefficients of 2 and 0.2 m 2 s −1 for momentum and tracers, respectively. Vertical mixing uses the K-profile parameterization (KPP) scheme of Large et al. (1994). The model is initialized from rest; the initial temperature and salinity distributions are uniform in the horizontal, and the vertical profiles are taken from the World Ocean Atlas 2013 climatology (Boyer et al., 2013). Initial fields of the biogeochemical tracers (O 2 , P and Z) are also horizontally homogeneous, and their vertical profiles are obtained from a water column (e.g. 1D) model version of the O 2 PZ model (not shown). The channel is periodic in the meridional boundaries; the eastern boundary (coast) is modelled as a free-slip wall, and the western one (open ocean) is open with radiative boundary conditions for momentum and tracers. Within the westernmost 50 km, a sponge layer is applied with an increased horizontal viscosity of 200 m 2 s −1 and a nudging timescale of 3 d in order to nudge all prognostic variables to their initial values over the offshore region. The bottom stress formulation is logarithmic with a roughness height of 0.01 m. The model is forced by a cyclic uniform wind field. To promote the destabilization of the upwelling front, we introduce a small spatial perturbation to the wind field, whose amplitude is less than 1 % of the total wind speed and which has characteristic wavelengths of 10 and 50 km in the alongshore and cross-shore directions, respectively.
In the coupled model, the time evolution of DO, P and Z concentrations are given by, respectively, where ( We distinguish here the terms "photosynthesis" from "carbon fixation" (Behrenfeld et al., 2008). Indeed, photosynthesis is the production of ATP (adenosine triphosphate) and NADPH (reduced form of NADP+; nicotinamide adenine dinucleotide phosphate) by the photosynthetic-electrontransport (PET) chain, and these two products are referred to as "photosynthate". By contrast, carbon fixation refers to the use of photosynthate by the Calvin cycle to produce simple organic carbon products independently of light. SP2015 chose to highlight these two stages in their model formulation by decoupling carbon fixation from photosynthesis. The respiration by phytoplankton u r is given by u r ( The constants ν and c 3 have the same functions as δ and c 2 in u r . A Monod function is a logical formulation to have respiration depending on oxygen concentration, since if the fluid environment becomes anoxic, then plankton cannot breathe so that respiration terms u r and v r vanish. The fourth term parametrizes bacterial respiration (oxygen loss due to natural depletion, e.g. remineralization). On top of the oxygen that is produced and consumed within marine ecosystems, a considerable part of dissolved oxygen eventually becomes gaseous and is then exchanged at the ocean-atmosphere interface, contributing to the estimations that more than one half of atmospheric oxygen is produced in the ocean (Harris, 1986). The remaining term of Eq. (2) is the air-sea flux of O 2 , which is the product of the gas transfer velocity V a , which depends on the square of the wind speed and on the Schmidt number (Wanninkhof, 1992), with the difference of the oceanic DO concentration and the solubility of O 2 in seawater.
In Eq.
(2), phytoplankton growth g([O 2 ], [P]) is given by a linear growth term α([O 2 ])[P] minus intraspecific competition γ [P] 2 , where γ is the intraspecific competition intensity. The per capita growth rate is a Monod function is a parametrization of nutrient availability based on the nitrate-density (ρ) relationships measured during an oceanographic cruise which surveyed the IPUS during upwelling season (the 2007 MOU-TON -MOdélisation d'Un Théatre d'Opérations Navalescampaign; see also Rossi et al., 2010Rossi et al., , 2013. This term is included here to circumvent the absence of nutrient limitation in the original O 2 PZ model of SP2015. We chose to add a parametrization instead of a new equation for the nutrients because this last option would necessarily increase the complexity of the model. SP2015 used a formulation where high oxygen increases phytoplankton growth but limits oxygen production by photosynthesis. The Warburg effect (Turner and Brittain, 1962) corresponds to the decrease of the photosynthesis rate at high oxygen concentrations so that the choice of f ( ] production makes sense. Moreover oxygen stimulates photorespiration which reduces photosynthesis yield. As soon as a water molecule (H 2 O) undergoes photolysis in the chloroplast of oxygenic photosynthetic organisms, O 2 and ATP (adenosine triphosphate) are produced through both photosystems I and II. The products of photosynthesis (ATP) then enter the Calvin cycle to convert carbon dioxide and other compounds into glucose, which is the food that autotrophs need to grow. As such, we can relate the growth term α ( . Regarding zooplankton (Eq. 3), its feeding efficiency is a function of oxy- , where η is the maximum feeding efficiency and c 4 is the half-saturation constant. Zooplankton mortality is given by µ [Z].
The parameters of the O 2 PZ model were adjusted so that its 0D version (without advection or diffusion; see Table A3) has a steady state at values representative of the IPUS (see Appendix A).

Simulations
The coupled model is initialized at rest. Biogeochemical tracers are initialized from vertical profiles (Fig. 3a) obtained by a water column 1D model of the O 2 PZ model that was run for a time long enough to reach stable vertical distributions (not shown). Using these profiles instead of climatological initial conditions allows for reducing the drift of the 3D coupled model (not shown). The initial temperature and salinity pro-   (Table 1) is designed to investigate the sensitivity of the coastal DO inventory to upwelling duration and community structure (diatom dominated or not). We choose the duration of upwelling-favourable wind instead of its magnitude as a control factor because it has been shown to be a better predictor of coastal hypoxia (Feng et al., 2012;Forrest et al., 2011;Zhang et al., 2018). Community structure is set through the growth regime of the primary producers (phytoplankton) with respect to nutrient input. Thus, a diatom-dominated regime is simulated using enhanced growth rates due to nutrient input, for the (E)nhanced runs, while a non-diatom-dominated regime is simulated with growth rates (N)eutral to and (L)imited by nutrient input. We test the sensitivity of our coupled system to the relative dominance of diatoms in the community, as they have been suggested to decline over the long term, being gradually replaced by dinoflagellates and cyanobacteria (Gregg et al., 2017).
The effects of upwelling season length are studied by running two simulations with four and nine upwelling favourable wind cycles, maintaining a total simulation time of 90 d and considering a diatom-dominated phytoplankton community. The alongshore wind profiles (Fig. 3c) are based on the cyclic wind pulses that are characteristic of the summer-early-fall upwelling favourable conditions in the IPUS, especially those observed during the MOUTON campaign: 10 d pulses with a maximum wind speed of 12 m s −1 (Rossi et al., 2013;Torres et al., 2003). The cross-shore wind component is zero in all cases. Despite the fact that "true" wind patterns observed over western Iberia are undoubtedly more variable (both in terms of direction and magnitude) than our idealized forcing inspired from a specific campaign, it provides the quickest development of a clear upwelling front generating small-scale instabilities. Note that sensitivity tests were performed (not shown) and revealed that the model responses are qualitatively similar among different wind strengths, given that the magnitude is high enough to ensure the development of the unstable upwelling front.
In situ observations of the IPUS, collected during the MOUTON 2007 campaign (Rossi et al., 2010(Rossi et al., , 2013 and compiled within the WOA13 (Boyer et al., 2013), are used not only to initialize but also to perform a comparison with the outputs of our simulations. The ECC simulation used the wind forcing and phytoplankton growth regimes that most resemble the IPUS region and therefore was chosen as the reference simulation.
We compared the full ranges of densities ρ, DO concentrations [O 2 ] and chlorophyll a concentrations [Chl a] of the reference simulation to all compiled measurements collected during the MOUTON campaign. Given the low complexity of the biogeochemical model and the highly idealized physical setting, a fully quantitative comparison is out of the scope of this paper; however we confirm that the coupled model reproduces qualitatively well the ρ-[O 2 ] and ρ-[Chl a] relationships obtained from in situ observations.
3 Results and discussion

Dissolved oxygen in the idealized IPUS
Our coupled model reproduces well the upwelling circulation and the typical biological responses. Indeed, the upwelling of nutrient-rich waters induced by the mean winddriven circulation promotes the growth of phytoplankton, leading to increased oxygen production by photosynthesis and the subsequent enrichment of oxygen in shelf waters (x < 80 km; Fig. 5a), as compared with the lower [O 2 ] found in offshore waters (x > 80 km). It also shows a low-[O 2 ] cell (∼ 50 mmol m −3 ), at the subsurface over the outer shelf (centred at 40 km and 60 m depth) because of the upwelling of low-O 2 waters, consistent with the shelf low-[O 2 ] cell (< 200 µmol kg −1 ) measured during the MOUTON cruise (e.g. Rossi et al., 2013). In a qualitative sense our idealized model results are similar to more realistic model studies, such as Gutknecht et al. (2013), which modelled the Benguela upwelling region and found a low-[O 2 ] plume for the climatological month of December at the shelf edge, and to observational studies such as Hales et al. (2006), which for the Oregon coast, measured [O 2 ] of 70-110 mmol m −3 in upwelled water at the shelf break, about 200 mmol m −3 less than at the surface, which is the same range of the vertical gradient simulated here.
The dynamics of the ECC simulation (Appendix B), are consistent with those documented by Durski and Allen (2005). A mean Ekman upwelling transport is established, moving subsurface waters upwards near the coast. The eddyinduced circulation (Cerovečki et al., 2009;Plumb and Ferrari, 2005) has the opposite sense to the mean circulation and is the strongest in the 50-to-100 km offshore range (Fig. B1). The result of the eddy-induced onshore and downward circulation is mainly seen in the subduction of the oxygen-rich waters just offshore the shelf edge, as suggested by the spatial pattern of the eddy stream function (Fig. B1).
The spatial coincidence between the high-[O 2 ] and highphytoplankton [P] concentration near the coast (Fig. 5)   toplankton growth is limited by the lack of nutrients in the euphotic zone.
The subsurface [P] maximum between 60 and 80 m depth results from the trade-off between sufficient oxygen, light levels and nutrients availability. Consequently, the enhanced phytoplankton growth rate due to nutrients appears to compensate for the low oxygen levels, thus promoting phytoplankton growth. Moreover, the low [O 2 ] will maintain the zooplankton stocks at low levels, thus constraining the grazing pressure of zooplankton on phytoplankton.
The dominant length scales, i.e. the inverse of the wave number where the alongshore [O 2 ] anomaly field spectrum has a maximum, vary with time and the distance to the coast (Fig. 6a). After day 30, between 50 and 150 km offshore, we observe the emergence of dominant length scales of 150 km in an intermittent fashion, in a process that agrees with previous studies of the growth and decay of large mesoscale structures in coastal upwelling systems (Durski and Allen, 2005).
Slopes of time-averaged spectra of alongshore [O 2 ] anomalies at three offshore locations (Fig. 6b) show, however, that the submesoscale has also a marked effect on the [O 2 ] field. All spectra are shallower than the k −3 slope of geostrophic turbulence (Charney, 1971), implying that the submesoscale is prominent in our simulations (Capet et al., 2008). The spectrum at x = 50 km is the shallowest and the closest to the k −5/3 slope in the range 8.10 −4 < k < 10 −4 . In the lower k range (k < 10 −5 ), the spectra at 50 and 100 km have higher energy than the spectrum at 200 km, as the dominant length-scale plot shows the largest dominant lengthscales in the 50-150 km cross-shore region.
The probability density function (PDF) of the [O 2 ] gradient norm (Fig. 6c) has a heavy tail distribution, associated with intermittency (Capet et al., 2008;Castaing et al., 1990). Turbulence redistributes oxygen through fronts and filaments (strain-dominated regions) as well as vortices (rotation-dominated regions). To analyse this process, we use the Okubo-Weiss criterion γ , which separates turbulent velocity fields in strain-dominated (γ > 0) and rotationdominated (γ < 0) regions. More specifically, we compute the daily average surface [O 2 ] for each γ and binned shore distance x.
For day 65, the [O 2 ] distribution by (γ , x) (Fig. 6d) Fig. 5a), in the range γ ∼ 0. Deviations from this background gradient are found along the cross-shore direction as we move to larger |γ |, implying that [O 2 ] anomalies with respect to the mean surface field are mainly found in filaments and vortices. Since the highest [O 2 ] levels are found in the shelf region, due to photosynthetic production, their presence further offshore strongly suggests offshore export by turbulent structures. We emphasize that offshore transport of [O 2 ]-rich water results in a net loss of O 2 for the continental shelf region.
The [O 2 ] anomalies found offshore in regions of positive and negative γ highlight the transport and dispersion of O 2 by (sub-)mesoscale structures from the shelf to the open ocean. We find that although the mean [O 2 ] as a function of γ (averaged over x in each γ bin) remains at ∼ 200 mmol m −3 (Fig. 6e), the [O 2 ] distribution is biased toward negative γ , as found by Combes et al. (2013) for tracers in the California current. Our results are also in line with those of Lovecchio et al. (2017) that demonstrate the importance of mesoscale processes for the zonal advection of tracers, with a key role of eddies, especially offshore the upwelling front.
The relative importance of the physical vs. the biological processes in controlling coastal oxygen dynamics was anal- ysed in terms of the volume-averaged rate of change of DO due to physics (advection and mixing), biology (photosynthesis, community respiration and remineralization) and airsea fluxes within a shelf control volume (from the coast to 50 km offshore and from the surface to 140 m depth) and averaged from day 20 to 90. The net enrichment rate, accounting for the physical and biological forcing, is positive (2.26 mmol m −3 d −1 ), and so, over time, the coastal oxygen inventory increases. DO concentrations are increased by biogeochemical processes at a rate of 4.27 mmol m −3 d −1 and by air-sea exchange (3.41 mmol m −3 d −1 ) but are limited by physical processes that remove O 2 from the coastal region (−5.42 mmol m −3 d −1 ) through lateral advection.

Sensitivity of dissolved oxygen to upwelling duration and phytoplankton growth regimes
The reference simulation (ECC) has shown three major features: the appearance of subsurface low-[O 2 ] cells because of the upwelling of O 2 -poor waters, the O 2 enrichment of surface levels, and the appearance of an [O 2 ] gradient between the biologically mediated enrichment at the inner shelf and the depleted offshore region. This gradient is, on the one hand, maintained by the upwelling circulation and, on the other hand, smoothed out by the cross-shore transport of O 2 by small-scale vortices and filaments. These effects are driven by the wind regime, which controls the upwelling dynamics and instabilities, and by the phytoplankton growth rate, which in turn affects phytoplankton concentrations. We now look at the sensitivity of the coupled model to these two drivers.
The continuous development of filaments and vortices in the ECC run maintains the offshore transport of dissolved oxygen (Fig. 7a). The wind shutdown in the ECS run causes a gradual decrease of turbulence which implies a diminution of the [O 2 ] redistribution by the cross-shore export physical processes, resulting in the creation of a sharper boundary between the oxygen-rich coastal waters and the depleted open ocean (offshore the upwelling front; Fig. 7b).
The phytoplankton growth regime impacts the mean [O 2 ] distribution, since nutrient neutrality (NCC) or limitation (LCC) results in the decrease of [O 2 ], when compared with the ECC run. The largest difference between the NCC and the ECC simulations (Fig. 8a) is observed in the coastal region, with maximum deviations from the ECC run of up to 100 mmol m −3 . Offshore the front, on the other hand, we find an increase of [O 2 ] with respect to the ECC run, with the difference reaching ∼ 50 mmol m −3 at x = 250 km. This positive difference is due to the absence of nutrient limitation so that phytoplankton growth and oxygen production persist in the nutrient-depleted surface levels.
The difference between the LCC and ECC cases (Fig. 8b) is also the highest in the shelf region, slightly above 100 mmol m −3 . Offshore the front (x > 50 km), the difference reduces to less than 50 mmol m −3 , above 80 m depth, vanishing below this depth. In both cases, there is a signature of the low-[O 2 ] upwelling over the slope and shelf edge and the eddy-induced subduction just offshore the front (x = 50 km).
When analysing O 2 budgets in the coastal control volume, from the shore to 50 km offshore and from the surface to 200 m of depth, the short upwelling duration simulation (ECS) results in a slight increase of the shelf [O 2 ] temporalmean net rate of change by 6 % with respect to the ECC simulation (2.26 to 2.40 mmol m −3 d −1 ; Fig. 9a). Although every term in the oxygen mean budget for ECS is smaller than for ECC, the physical-sink decrease is slightly higher than the combined decrease in the biological source and sink (Q O 2 ) and air-sea exchange.
The changes in the physical sink also occur faster than those in the biological terms. The wind cessation at day 40 causes an immediate decrease in the physical term (Fig. 9b, lower panel), while only after day 70 can we observe a significant difference between the biological source and sink terms of both the ECC and ECS simulations.
Changes in the phytoplankton growth regime modify the physical-biological balance of oxygen in the coastal control volume. The coastal region is net autotrophic (Q O 2 > Overall, the phytoplankton growth regime has a stronger impact than the wind regime in the DO inventory in the shelf control volume. The NCC and LCC cases result in net rates of change of dissolved oxygen that are 75 %-80 % lower than in the reference case.

General discussion and conclusions
We built a coupled physical-biogeochemical model of an idealized seasonal coastal upwelling to study the effect of the upwelling season length and phytoplankton community structure on dissolved-oxygen inventory. The lowcomplexity O 2 PZ model accounts for oxygen production by photosynthesis and consumption by both respiration and remineralization, but its simplicity does not allow for including other oceanographic processes that could influence oxygen cycling over shallow shelves, such as for instance the complex benthic-pelagic coupled processes (Capet et al., 2016). While we used data from the Iberian Peninsula upwelling to initialize and compare the range of values simulated by our model, our idealized configuration allows for drawing general conclusions about the mechanisms governing the dissolved-oxygen levels over the continental shelf. When compared to measurements, our model reproduces qualitatively the O 2 -density relationship as well as the upwelling of oxygen-poor waters onto the shelf and the offshore transport of oxygen due to filaments and vortices. While the addition of air-sea exchange processes as well as our novel nutrient-density parametrization have probably contributed largely to simulate such realistic outputs despite the simple biogeochemical formulation, the coupled system could be further improved. One probably important task, which we keep for future studies, is to include realistic subsurface O 2 inputs from oceanic currents, as it has been shown to control largely O 2 dynamics over the shelf (Montes et al., 2014;Reboreda et al., 2015).
Although there is a clear separation of coastal waters with high [O 2 ] and offshore waters with low [O 2 ], the turbulent lateral transport of high-[O 2 ] water disrupts this picture, as oxygen-rich upwelling filaments and vortices are generated at the front and then travel offshore. This creates a sink of dissolved oxygen on the shelf that must be compensated by photosynthetic production and air-sea fluxes. In our reference simulation, with a repeated upwelling wind cycle, the net [O 2 ] rate of change was positive, owing to the predominance of the photosynthetic oxygen production and airsea fluxes over the lateral transport of oxygen, which acts to remove it from the shelf region. Considering a shorter upwelling season (the ECS case with four 10 d upwelling wind cycles followed by wind cessation) reinforces somewhat the oxygen enrichment trend. When the phytoplankton community is dominated by groups that show nutrient limitation or neutral growth (LCC and NCC cases), the biological source weakens, and the shelf becomes net heterotrophic. The physical sink is also affected, producing a weaker loss of dissolved oxygen in the coastal upwelling. Air-sea exchange then maintains the (much lower) net oxygen enrichment trend.
Our findings that the oxygen inventory net rate of change is inversely proportional to the duration of the upwelling season is consistent with other studies (Adams et al., 2013;Siedlecki et al., 2015;Zhang et al., 2018). In the present, the IPUS is not hypoxic, but, as the upwelling season duration in the IPUS is thought to increase in the future (Miranda et al., 2013;Wang et al., 2015), it suggests an increased vulnerability of coastal regions to an episodic or long-term decrease of DO levels. The result of lower phytoplankton growth rates driving the shelf region from autotrophy to heterotrophy also points to an increased risk of loss of dissolved oxygen if trends for community shifts towards smaller or less-efficient phytoplankton (Gomes et al., 2014;Gregg et al., 2017;Head and Pepin, 2010;Zhai et al., 2013) are specifically confirmed for the IPUS. In permanent oxygen minimum zones (OMZs), changes in physical forcing (e.g. wind regime) and/or biological characteristics (e.g. size structure or functional types) can have dire implications for marine ecosystems and fishery management. Indeed, the measured expansion of the tropical Atlantic and Pacific OMZs (Stramma et al., 2008) implies a shallower upper boundary of the low-O 2 core and an increased sensitivity of oxygen levels over the continental shelf to upwelling winds, which are indeed predicted to intensify with climate change (Sydeman et al., 2014). Future work could be focused on understanding this high sensitivity and on anticipating how it would affect the stability of the coupled system in a changing climate.

Appendix A: Determination of O 2 PZ model parameters
The biogeochemical component of the 3D coupled model (Eqs. 1-4) is a system of ordinary differential equations (ODEs) that computes the local time evolution of [O 2 ], [P] and [Z]. The ODE system accounts for the processes that move O 2 , P and Z between different reservoirs. Each ODE of the model then expresses the input and output of the species from the specific reservoir. The dynamical properties of the ODEs will therefore influence the global behaviour of the coupled physical-biogeochemical model. When using coupled models such as Eq. (1) (2)-(4) that does not change with time and attracts nearby system trajectories. This guarantees that, when removing the advection and diffusion terms, the system will evolve to a constant state whose properties are already established (e.g. positive concentrations). The fixed points of the O 2 PZ system are the states X + , where (A1) Since the system is nonlinear, a solution for Eq. (A1) is not straightforward, and approximate methods are then a good choice. Since the model will be used to study O 2 dynamics in coastal upwelling systems in the Iberian Peninsula western shelf and slope, it is logical that the model uses initial conditions that are characteristic of this region. The characteristic state X w = ([O 2 ] w , [P] w , [Z] w ) may be taken from a climatological dataset, in situ observations, remotely sensed data, etc. as long as it is representative of the actual values of ([O 2 ],[P], [Z]) in that region. Regarding the position of X w in the state space, a simple proposition is that X w is close to a coexistence steady state X 0 . This means that, if the model is initialized at X w , it will converge to X 0 and will not evolve to an extinction state or other steady states with unrealistic or uncharacteristic values of ([O 2 ],[P],[Z]) for the region. This is the approach adopted here to determine the values of the parameters of the system, and we adopt the characteristic values shown in Table A1 for the Iberian Peninsula upwelling.
The concentration of DO is taken from the climatological mean at 12 • W, 41 • N and 10 m depth. The phytoplankton concentration is obtained from ocean colour images of Chl a in the western Iberian shelf at approximately the same position of the DO reading, after converting it to nitrogen concentration using a ratio of 1.59 between Chl a and nitrogen content (Montes et al., 2014). Finally, the zooplankton concentration is found by using a phytoplankton-to-zooplankton ratio of 3.6, taken from SP2015 for the case of a stable coexistence state of the nondimensional model. For some parameters we adopt the values that are used elsewhere in the literature. To find the most adequate values of the remaining parameters, we search for the parameters that make possible the propositions above regarding the stable nature of X w . There are 11 parameters that relate O 2 with phytoplankton and zooplankton. It is a considerable amount of free parameters to handle. To constrain the choice of values for this set of parameters, making it easier to tune the model, we introduce the following three assumptions, regarding the model behaviour at X w : 1. At X w , zooplankton growth should not be hampered by the value of [O 2 ] w . This requires that the feeding efficiency term is at or near its maximum h.
2. At X w , the phytoplankton growth rate should be near or at its maximum value of B.
3. At X w , phytoplankton O 2 production rate is significantly larger than O 2 respiration rate.
These assumptions are adopted to facilitate a coexistence state at X w , by not limiting the ability of zooplankton to feed on phytoplankton and at the same time by facilitating phytoplankton growth to feed zooplankton and to produce O 2 in a greater amount than it respires. By varying the A, B and h parameters, the fixed point of the system can be made to approach X w . Although other parameters may be varied instead of these three, we found that these are the ones whose variations are easier to judge in terms of displacement of the isoclines. This results in a coexistence fixed point X 0 that is located near X w , with coordinates shown in Table A2, which is a stable spiral point with eigenvalues of the Jacobian λ 1 = −0.0847 and λ 2,3 = −0.0009 ± 0.0221i. We analyse the dynamics of the ECC simulation in terms of turbulent energy transfers, mean and eddy-induced circulations, and dominant turbulent length scales. Mean quantities refer to the alongshore average, e.g. for the cross-shore velocity u: u(x, z, t) = L −1 u(x, y, z, t)dy, where L is the alongshore length of the channel. Perturbations are the departures from this mean: u = u − u and u = 0. The mean kinetic-energy transfer term is where ρ 0 is the reference density, v and w are the alongshore and vertical velocities, and a subscript denotes partial differentiation. The perturbation potential energy transfer term is c pe = −g w ρ .
In the ECC case, c pe remains positive, increasing after the wind intensification part of the wind cycle and decreasing when the wind relaxes (Fig. B1a), while c ke oscillates between positive and negative values. The mechanism for c ke growth seems to be the coalescence of disturbances into a single large eddy that is subsequently perturbed by rapidly emerging small-scale patterns (Durski and Allen, 2005). Anti-correlation between c pe and c ke is observable at days 35, 47, 53 and 60 (Fig. B1a), suggesting that c ke evolution is related to nonlinear wave-wave interaction rather than to wave-mean flow interaction (Durski and Allen, 2005). For sustained winds, like the ECC wind cycle with short wind intensification and relaxation periods, the effects of this interaction were significant since the start of the simulations (Durski et al., 2007).
Initially, the surface density peak length scale l e , i.e. the most energetic scale, is homogeneous, and l e is equal to the domain length (Fig. B1b).
Shortly after the start, small-scale (O(10 km)) fluctuations appear close to the shore, then grow in length and spread offshore. From day 20 onwards, l e fully develops and ranges from 20 to 160 km across the shore, although until day 36 there are still regions where density is rather homogeneous.
Offshore the upwelling front (roughly from 50 km) we observe the emergence of two distinct regions in the crossshore direction: a region adjacent to the front where largescale structures are often dominant (l e ∼ 150-160 km) and further offshore where smaller-scale structures are dominant (l e < 100 km). The cross-shore extent of the first region appears to be related to the excursions in the c ke time series (Fig. 3a), judging by the coincidence of the largest fluctuations of c ke with the widest swaths of l e ∼ 150-160 km (days 45, 60, 68 and 80). The second region further offshore seldom exhibits l e > 100 km, as it is usually populated by filaments emerging from the front with contrasted density and by vortices resulting from the roll up of filament tips. The wind-induced mean circulation in the x-z plane is = −τ/ρ 0 f , where τ is the applied wind stress and f is the Coriolis parameter. The eddy stream function is ψ = ( v b b z −λ 2 w b b y )/(λ 2 w b b y + b z ), where b is buoyancy and λ = 1000 is a stretching constant that accounts for the large aspect ratio of oceanic flows (Nagai et al., 2015;Plumb and Ferrari, 2005). The wind-induced mean circulation in the vertical cross-shore plane (Fig. B1c) is stronger at the surface and near the coast and weaker in the interior. A typical upwelling circulation is simulated: Ekman transport moves surface waters offshore and upwards near the coast, giving rise to the mean upward tilting of the isopycnals. The eddy-induced stream function (Fig. B1d) represents the component of the flow across the mean isopycnals (Cerovečki et al., 2009;Plumb and Ferrari, 2005), opposes the mean circulation and can surpass it. In agreement with Cerovečki et al. (2009) we find that the eddy stream function is stronger at the diabatic surface layer, in a thin layer (up to 20 m depth) with strong shoreward circulation 100 km offshore, deepening to 150 m at 75 km offshore. There is another strong eddy downwelling cell, on the shelf-slope transition, that goes down to 160 m depth.
The submesoscale and mesoscale turbulence is sustained by energy transfer from the potential energy field to the turbulent kinetic-energy field, through baroclinic conversion (Durski and Allen, 2005;Marchesiello et al., 2003). In agreement with previous modelling studies (Capet et al., 2008;Durski and Allen, 2005), the global picture is thus that the intermittent wind bursts act as a recharge for the mean potential energy which is then released to eddy kinetic energy through baroclinic instability. When the wind decreases, the eddy flux induces a rapid re-stratification of the upper layer (Capet et al., 2008).
Code and data availability. Simulated data and codes are available upon request.
Author contributions. JHB, VR and VG designed the study with input from LR, YM and PH. JHB performed the simulations. JHB, VR and VG analysed the results and prepared the first draft of the paper. All authors critically reviewed, commented and prepared the final paper.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. This work is part of the TEASAO project, which is funded by the "IDEX Attractivity Chair" programme of the University of Toulouse. João H. Bettencourt acknowledges postdoctoral fellowship financial support from the TEASAO project and CENTEC (University of Lisbon) for hosting him during his stays. Vincent Rossi acknowledges fruitful discussions with Martinho Marta-Almeida and travel support from INSU-CNRS. The MOUTON 2007 campaign was carried out as part of a SHOM (Service Hydrographique et Océanographique de la Marine) project, and we acknowledge the competence and investment of the SHOM technical staff who took part in it.
Financial support. This research has been supported by the IDEX UNITI -University of Toulouse (TEASAO IDEX UNITI -University of Toulouse).
Review statement. This paper was edited by Ana M. Mancho and reviewed by two anonymous referees.