High-resolution modelling of long-term trends in the oxygen and carbon cycles of the Benguela upwelling system

We investigate driving forces of the biogeochemistry of the Benguela upwelling system (BUS) and their temporal changes over the 20th century. For this purpose, we developed a global ocean-only model in a stretched grid configuration, which resolves meso-scale circulation structures in the area of interest. The biogeochemical module of this model is extended by a more comprehensive nitrogen cycle to account for the specific nitrogen loss processes common in eastern upwelling systems. The model is forced by 110 years of atmospheric reanalysis data. To assess the impact of meso-scale circulation structures 5 on local biogeochemical processes we compare our results to a set-up with a coarser horizontal resolution, comparable to the ocean component of Earth system models used for anthropogenic climate projections. In the higher spatial resolution we find enhanced intermediate depth ventilation (200-1000 m) with concurrent reduced loss of bioavailable nitrogen and a high shelfbound biological production, in line with observations. Moreover, only in the high resolution setup do multi-decadal trends of deoxygenation match observation-based estimates. Our study supports the view that the presence of meso-scale circulation 10 structures exerts a major influence on biogeochemical patterns, especially on mid-depth oxygen concentrations. Furthermore, we show for the first time that by including this high spatio-temporal variability of the circulation, the regional anthropogenic carbon uptake of the BUS over the 20th century is lower than in the coarse resolution model. This indicates that, at least for some regions, the pathway of changes in the marine biogeochemistry as projected by state-of-the-art coarse resolution Earth system models is associated with high uncertainty. 15 1 https://doi.org/10.5194/bg-2022-27 Preprint. Discussion started: 10 February 2022 c © Author(s) 2022. CC BY 4.0 License.


Introduction
There is increasing interest in the role of coastal upwelling regions within the climate system of the Earth (Gruber et al., 2012;Emeis et al., 2018). These regions, and especially the Benguela upwelling system (BUS), host the most productive ecosystems with corresponding high organic carbon fluxes to the ocean interior. Decay of organic material in combination 20 with the hydrodynamical characteristic of a low ventilation of intermediate depths (200 to 1000 m) creates the world's largest oxygen minimum zones (OMZs) (Cabré et al., 2015). Furthermore, high organic fluxes, often referred to as biological carbon pump, also contribute to the global marine carbon storage. Climate induced changes in coastal upwelling systems put multiple stresses on their ecosystems which might cause unprecedented impacts on the global carbon cycle (Gruber et al., 2012;Bopp et al., 2013). It is hypothesised that ocean warming increases deoxygenation (Schmidtko et al., 2017) which could affect nutrient 25 cycling and the biological pump with potentially adverse consequences on marine ecosystems. Increasing levels of atmospheric CO 2 lead to a lowering of seawater pH. Having already naturally low levels of seawater pH, the biological ecosystem in coastal upwelling regions are particularly vulnerable to ocean acidification (Gruber et al., 2012). Thus, identifying controls of biogeochemical processes in coastal upwelling systems is essential to assess the potential impacts on deoxygenation and regional carbon budgets in a changing climate.

30
A common feature of eastern boundary upwelling systems, including the BUS, is their high spatial and temporal variability in physical processes in response to coastal jets and a quasi-permanent wind-induced upwelling modulated by a wind stress curl driven current system (Fennel, 1998). Assessing mean state and temporal evolution of the biogeochemistry in coastal upwelling systems is, therefore, an ambitious task. The high spatio-temporal variability of the physical drivers challenges observational methods to achieve a comprehensive picture. Modelling approaches could give insights in the physical and biogeochemical 35 pattern of the BUS. However, none of the global Earth System models, which are currently used for future climate projections, is capable to reproduce present-day oxygen concentrations in the BUS (Séférian et al., 2020) or past decadal trends (Buchanan and Tagliabue, 2021). They overestimate the size of the OMZ in the eastern South Atlantic near the coast of South Africa. The simulated multi-model mean of the decadal O 2 trend is far too weak because models show no consistency in the trend (see Fig. S1 in Buchanan and Tagliabue (2021)). Séférian et al. (2020) postulated a systematic bias in the ocean biogeochemical 40 models which seems to be independent from ocean resolution or complexity of marine biogeochemical models. However, the horizontal nominal resolution of the considered Earth system models is typically around 100 km; only one could be classified as eddy-permitting (25 km). Therefore, it is debatable if the simulated bias in the OMZs is not attributable to the low spatial model resolutions which do not resolve meso-scale structures. These structures, being eddies and filaments, are known to introduce strong ocean mixing and water column ventilation. Recently, a global eddy resolving coupled ocean model (spatial Philander (1981) and a wind-driven turbulent mixing scheme for the upper ocean. Details on parametrizations of atmospheric exchange fluxes can be found in Jungclaus et al. (2013). 85 Here we use MPIOM in two different grid configurations, both covering the entire world ocean. GR15 is the standard ocean model component used within the MPI Earth System Model (Jungclaus et al., 2013;Mauritsen et al., 2019). GR15 is a bi-polar grid with one pole is located on Greenland and one on Antactica. It has a horizontal resolution of 1.5 degrees gradually varying between 15 km near Greenland and more than 180 km in the tropical ocean. The second setup uses a special grid configuration (higBUS, Fig. 1) where the poles are shifted to South Africa (27 • S, 24 • E) and South America (20 • S, 58 • W). Shifting of 90 the pole locations leads to a finer horizontal resolution between 10 • S-40 • S in the Atlantic (see e.g. Mathis et al. (2019) for the same method applied to the North Atlantic). The Atlantic sector of the Southern Ocean (>60 • S) has a similar horizontal resolution in both setups. The location of the poles leads to a fine grid resolution in the area of interest but a relatively coarse resolution of approx. 330 km in the Pacific. The Benguela current system is resolved with approx. 10 km, a resolution which can be considered as eddy resolving. The Rossby radius of deformation for the BUS is estimated to be about 30 km (Hallberg,95 2013; Chelton et al., 1998). In MPIOM, the impact of the parameterization for unresolved eddy-induced transport by Gent et al. (1995) is scaled with the horizontal grid distance, basically turning off this artificial eddy transport in regions such as the BUS in higBUS. The vertical resolution is identical in both setups with 40 layers. In the upper 100 m, layer thicknesses are 4 https://doi.org/10.5194/bg-2022-27 Preprint. Discussion started: 10 February 2022 c Author(s) 2022. CC BY 4.0 License. around 10 m and they are gradually increasing towards ocean interior (layer thickness max. 600 m). The topography is based on ETOPO2 (NOAA, 2006) and we consider partial grid cells in the last wet grid cell. The time step is 1 hour in GR15 and 30 100 min in higBUS. Tidal effects are not included.

Biogeochemical model
One focus of the global biogeochemical model HAMOCC is and has always been on a good representation of large scale geochemical tracer distributions (Maier-Reimer, 1993). Aspects of the marine biota have been kept as simple as possible. This version of HAMOCC considers effects of the marine biota by an extended plankton dynamics module representing nutrients 105 (phosphate, nitrogen species, iron), two types of phytoplankton (bulk and diazotrophs) and one type of zooplankton, detritus and dissolved organic matter (DOM) (Six et al., 1996;Ilyina et al., 2013;Paulsen et al., 2017). Shell material production in forms of opal or calcium carbonate is linked to the production of detrius with a preference to opal production if silicate is available. Organic matter/shell material production and decomposition/dissolution affect the alkalinity. All biogeochemical tracers are subject to the same hydrodynamics as temperature and salinity. Particulate matter such as detritus and shell material 110 are additionally subject to gravitational sinking. As long as oxygen is available detritus and DOM are remineralized which releases carbon and nutrients in a fixed composition ratio ("Redfield ratio", Takahashi et al. (1985)) and consumes oxygen and alkalinity. In suboxic zones denitrification and sulfate reduction control the decay of detritus. Particles reaching the sea floor fuel a sediment module with 12 active layers (Heinze et al., 1999). Each sediment layer consists of solid particles and porewater which interact following the same processes of remineralization, denitrification and sulfate reduction as represented 115 in the water column. Diffusive exchanges between the sediment layers and the water-sediment boundary supply oxygen and nitrate to fuel organic matter decay and maintain the reflux of carbon, phosphate, iron, silicate and alkalinity from the sediment to the water column. Depending on the particle flux and the sediment composition, shells and organic material might be buried and transferred to deeper sediment layers and eventually end up in the diagenetically consolidated burial layer.
Particle fluxes of organic matter are described by the linearly increasing sinking speed following the formulation of Martin 120 et al. (1987). In the upper 90 m the sinking speed is kept constant at 3.5 m d −1 and it reaches a maximum of 70 m d −1 at 6000 m. Shell material is sinking with a constant speed of 30 m d −1 for calcareous shells and 22 m d −1 for opal frustules. Processes of shell dissolution and air-sea gas exchange of CO 2 are formulated as in Ilyina et al. (2013).
A new development that we are introducing for the HAMOCC version used here is a more comprehensive nitrogen cycle, which in particular includes a detailed representation of nitrogen loss processes. Nitrate reduction processes occur only in 125 suboxic to anoxic waters. They potentially lead to the loss of bioavailable nitrogen and, thus, affect directly the biological production. In the following we concentrate on the description of newly implemented processes of the N-cycle and give only a general overview of other represented biogeochemical processes.

Nutrient cycle
In general, biological activity consumes nutrients in the euphotic surface layer to build organic material. This is partly rem-  (2)-(15). and the sediments. The local change in nutrient concentrations occurs due to hydrodynamics, phytoplankton growth, grazing activities, organic matter decay, fluxes to the sediments and weathering input. In case of total inorganic nitrogen (TIN) which comprises nitrate, nitrite, and ammonium we also have to consider air-sea gas exchange for ammonia. In general, all biogenic tracers follow where U is the three-dimensional advection vector, D the diffusion operator, and Ψ i k (C k , C i ) refers to the source and sink operators for biogeochemical processes, fluxes across the air-sea interface (Ψ airsea ) and interactions with the sediment (Ψ sedf lux ).
In the following paragraphs we describe the processes summarized by Ψ in equation (1) on TIN. A schematic view on the processes is given in Fig. 2. All parameters and rates used in the new N-cycle of the model are given in Table A1 in the Appendix.

140
A detailed description for other prognostic tracers like phosphate, dissolved iron, dissolved inorganic carbon, and alkalinity is found in Ilyina et al. (2013).
Sinks for NO − 3 are phytoplankton growth, and dissimilatory nitrate reduction to nitrite (DN RN ) or ammonium (DN RA). Sources for NO − 3 are oxidation of NO − 2 (O NO 2 ) and as a small byproduct of anammox (Anam). Correspondingly, changes of NO − 2 and NH + 4 are given by where N RN 2, the denitrification of NO − 2 to N 2 , is a total N-loss process of bioavailable nitrogen. G phy is the nutrient consumption due to phytoplankton growth on NO − 3 or NH + 4 . The NH + 4 source from remineralization includes decomposition of organic material via remineralization in oxygenated waters (F det ), grazing activity on phytoplankton F pz , grazing activity on zooplankton M zoo , and bacterial decomposition of dissolve organic matter λ dom DOM. In HAMOCC, we use a constant ele-165 mentary composition of organic material following the "Redfield" ratio R N :P of P:N:C:-O 2 = 1:16:122:172 (Takahashi et al., 1985). Phytoplankton production is limited by the availability of nutrients (phosphate, nitrate, ammonium or iron), light and temperature. For a detailed description of the standard dependencies of growth on light and temperature we refer to Ilyina et al. (2013).
For each nutrient, we determine a potential phytoplankton growth G phy with the corresponding half saturation constant for that 170 nutrient separately and, then, use the minimum 7 https://doi.org/10.5194/bg-2022-27 Preprint. Discussion started: 10 February 2022 c Author(s) 2022. CC BY 4.0 License.
for i being P or Fe Phytoplankton favours the uptake of ammonium, if present, which is considered by the inhibition factor (Fennel et al., 2006) σ inhib = 1 1 + N NH 4 /bk NH 4 J(T, I(z)) is the maximal growth rate. It is depending on temperature (Eppley, 1972) and available photosynthetic radiation 180 at the depth z (see Paulsen et al. (2017) for more details).
The process descriptions for grazing activities and bacterial decomposition are identical to Ilyina et al. (2013), but aerobic and anaerobic remineralization processes are updated. Aerobic remineralization of F det is a linear function of detritus concentration occurring only in water with an oxygen concentration above 2 mmol m −3 .
All nitrogen in organic material is remineralized to ammonium and, thereby, less oxygen is consumed than it is needed to remineralize organic matter to nitrate (Paulmier et al., 2009). In oxygenated water stepwise nitrification converts NH + 4 and NO − 2 back to NO − 3 . These processes are inhibited in the presence of light (Fennel et al., 2006).
with σ I = bk I (bk I + I(z)) In oxygen minimum zones (OMZs) , i.e.suboxic waters with oxygen concentrations below O 2crit , processes of dissimilatory nitrate reduction to nitrite (DN RN ) or ammonium (DN RA), and heterotrophic denitrification (N RN 2) degrade organic  in the presence of oxygen. Only N RN 2 leads to a complete loss of fixed nitrogen while organic material is degraded.
The second pathway for complete loss of fixed nitrogen is the anaerobic oxidation of NH + 4 with NO − 2 , anammox. Anammox bacteria combine ammonium and nitrite to dinitrogen and hereby produce small amounts of nitrate (Brunner et al., 2013).
The concentration of NH + 4 is generally low in the ocean which means that the magnitude of Anam is mainly controlled by 205 the availability of NH + 4 . σ ano2 is defined to be equal 1 at an oxygen concentration larger than the threshhold value O 2crit and it decreases with decreasing oxygen according to This allows that DN RN, DN RA, and Anam start to occur in suboxic waters as it has been observed .
Complete denitrification N RN 2 is restricted to oxygen concentrations below 2 mmol m −3 . All rates λ and half saturation 210 constants bk are given in Table A1.
Please note, that all processes described above are active in the water column and the sediment. Diffusive fluxes for all nitrogen components at the water-sediment boundary are also included.

Air-sea gasexchange of ammonia 215
Ammonia (NH 3 ) is a soluble gas and its air-sea transfer is generally considered to be under gas-phase control (Liss, 1983).
The concentration of NH 3 is determined by the concentration of NH + 4 , the dissociation coefficient for ammonium, which is a function of temperature (T ), salinitiy (S), and the ionic strength of sea water H + .
with 220 K a = 10 −pKa pK a = −0.467 + (0.00113S) + (2887.9/T ) The constants in Eq.(14) are taken from Johnson et al. (2008). The flux of ammonia across the ocean-atmosphere interface is a function of the gas-phase transfer velocity k g and the concentration gradient for NH 3 between ocean and atmosphere. Here, we assume that the atmospheric concentration of NH 3 is neglectable (atmospheric concentrations are generally less than 1 nmol m −3 , given no terrestrial sources) and the flux is only a function of the sea water concentration times the dimensionless coefficient of Henry's law for ammonium K H .

Nitrogen fixation
The loss of fixed nitrogen by processes of denitrification and anammox in the OMZ is compensated by the fixation of elemen- Tg N yr −1 ) for present day conditions. Here, we are briefly summarizing the major features of this prognostic treatment and we refer to the full description of this model development to Paulsen et al. (2017).
As for bulk phytoplankton, diazopthrophs live in the upper ocean controlled by the avalability of nutrients, light and temperature. In contrast to bulk phytoplankton, diazothrophs, especially the here considered Trichodesmium, have a distinct optimal temperature range in which they can grow. This limits their distribution to the zonal band of 40°N to 40°S (Paulsen et al., 2017).

245
Diazotrophs are described to have a smaller maximal growth rate than bulk phytoplankton, but they have the same biochemical composition of their organic tissue following the fixed Redfield ratio. In regions with sufficient nutrient supply diazotrophs have to compete with bulk phytoplankton. Their ability to fix N 2 for their N supply allows them to also grow in regions which are already depleted with respect to bioavailable nitrogen compounds. The mortality of diazotrophs is considered a first-order process with a constant decay rate. The major part of the organic material is directed to the detritus pool (90 %) and subject to 250 export by gravitational sinking. The remaining part is consider as dissolved organic matter which decays with a constant rate (Paulsen et al., 2017;Ilyina et al., 2013).

Model setup and transient simulations
Our goal is to run a simulation with a consistent forcing over the 20th century allowing us to analyse decadal trends. To this 255 end, we need to perform a spin-up run which, in its final state, has a low artificial model drift and shows a reduced initial shock when applying the transient ERA20C forcing.
To achieve starting points for the transient simulation for both model setups, a spin-up procedure was done in several steps.
First, we forced both models with the same atmospheric boundary conditions from a climatology based on the second ECMWF Re-Analysis project (ERA-40) (Simmons and Gibson, 2000) from which a mean annual cycle in daily resolution for relevant In a second step, we switched to the ERA20C forcing, which is ECMWF's first atmospheric reanalysis of the 20th century, from 1901-2010 (Poli et al., 2016). We converted the data from a 3-hourly to a 6-hourly forcing. The continental runoff was adopted from the climatological OMIP forcing (Röske, 2006). For the spin-up phase we use only the years 1910 to 1939 of the ERA20C forcing as a cyclic forcing. We prescribe a fixed atmospheric pCO 2 corresponding to the year 1900 (295.7 ppmV).

280
After about 960 years we started the transient simulation from 1901 to 2009 with the corresponding ERA20C forcing and increasing atmospheric pCO 2 following the historical reconstruction (https://data.giss.nasa.gov/modelforce/ ghgases).
We are aware of the fact that a good spatial resolution of the local wind stress field is essential to capture some of the main  1979 -present) and therefore includes already a strong anthropogenic signal. Other reanalysis products, such as NOAA-CIRES-DOE 20th Century Reanalysis version 2 (http://www.esrl.noaa.gov/psd/data/20thC_

Additional Diagnostics
To gain more insight into the flow pattern and water mass composition we introduce artificial dye tracers. Dyes are passive tracers and follow the identical hydrodynamical processes as biogeochemical tracers. They have no sinks or sources within the ocean interior. Their boundary condition is set to 1 at the surface within their source region and set to 0 at the surface outside 295 of this region.
with ∂D i ∂dt = 1 at the surface within source region of D i We define 9 source regions ( Fig. 3) to cover the whole model domain (without any overlaps). The relative contribution of a dye 300 concentration to the sum of all dye concentrations within a grid box tells us the regional origin of the water mass. In addition, we include a linear age tracer. It is set to zero at the ocean surface and treated as the dye tracers while constantly ageing in the ocean interior.
Furthermore, we track the advection of oxygen, phosphate, and detritus to investigate the eddy contribution versus mean flow transport. In general, the eddy transport for a tracer B at each geographical position can be estimated from the relationship U corresponds to the flow directions in zonal, meridional and vertical direction (u,v,w respectively). By substracting the temporal average of the product of the monthly mean velocity U times the temporal average of the monthly mean tracer concentration B from the temporal average of the monthly mean product U B we obtain the nonlinear part of the flow field which can be interpreted as eddy-induced transport (von Storch et al., 2016), but which includes interannual variability as well. 310 We use the term "eddy-induced" for U ′ B ′ even in the coarse resolution GR15, where meso-scale structures are not resolved.
Our analysis of meso-scale structures is based on monthly mean values because of the strong seasonality in the BUS.

General circulation pattern
The large scale circulation pattern and long term mean transports through major passages are similar between higBUS and GR15 and compare reasonably well to observations (Table 1). Due to the relative coarse resolution of the northern part of the 320 Atlantic Ocean and the Arctic the Atlantic Meridional Overturning Circulation (AMOC) in higBUS is rather low.
For the water mass composition of the BUS two prominent circulation patterns are relevant. The southern part of the BUS is strongly influence by the Agulhas current and its retroflection. Warm and salty water is spreading around the southern tip of Africa, turns southward and flows back into the Indian Ocean as indicated by the barotropic streamfunction of higBUS 325 (Fig. 4a). Part of this water, transported from the Indian Ocean northward into the Atlantic Ocean, feeds the Bengula current system (Agulhas leakage, Biastoch et al. (2009)). It has been already shown that eddy permitting model resolutions give a better representation of the Agulhas retroflection (e.g. Jungclaus et al. (2013)). The barotropic streamfunction shows two distinct cells  in the zonal band of 30-45 • S in higBUS being in line with observation based estimates (Wunsch, 2011). In contrast, GR15 shows only one "super cell". Models with coarser resolution tend to have a more laminar flow field with a higher transport of 330 Indian Ocean water into the Atlantic and a stronger Agulhas leakage at the tip of Africa (e.g. Marcello et al. (2018), Fig.1c ). This is underlined by the distribution of the dye tracer from SIO at 50-150 m in the South Atlantic (Fig. 5).
In higBUS, SIO water reaches the South Atlantic and is carried northward with the Benguela current without much westward spreading. Water from the SIO is the dominant fraction of the water mass composition south of 30 • S at 50-150 m in our simulation. Here, at the southern tip of Africa, water from the Agulhas leakage and South Atlantic water mix to form Eastern 335 South Atlantic Central Water (ESACW, Mohrholz et al. (2008)) as concluded from observations. The simulated dye tracer This water from SACW enters the BUS region at its northern boundary (Mohrholz et al., 2008) and is advected by a strong 340 poleward coastal current, very prominent at 23 • S (Fig. 6c), which pushes the branch of the Benguela current westward and prevents the spreading of SIO water along the coastal shelf (Fig. 5a, and Fig. A2). In higBUS, the poleward current is more intense and narrower than in GR15 with climatological mean velocities of more than 6 cm/s in the core at 100-200 m depth.
Close to the coast, we find a northward current which was also identified from current measurements at a mooring station at 23 • S and 14 • E (see Mohrholz et al. (2008), Fig. 14). Fig. 6a,b already give some indications of a more realistic cross shelf 345 exchange in higBUS where we find a zonal mean velocitiy of > 2 cm s −1 at the shelf break.
A water mass analysis by Mohrholz et al. (2008) for the same location (mooring data in 120 m at 23 • S, 14 • E) also supported a dominant fraction of SACW compared to ESACW (i.e. equivalent to the sum of SIO and SATL in our setup, see Fig. 3).
The study by Tim et al. (2018) tracking Lagrangian particles found a similar distribution with a high contribution of SIO water in the southern BUS (between 27 • S to 34 • S) and a lower presence of SIO in the northern part (between 17 • S and 27 • S).

350
Results of higBUS are consistent with these estimates, while GR15 is not able to resolve these complex circulation structures showing a different pattern with higher SIO and SATL fractions in the northern BUS (Fig. 5,Fig. A2).    To summarize, the higBUS simulation shows that the incorporation of meso-scale activity enables a reasonable representation of the major hydrodynamical pattern of the BUS. Discrepancies as described above still exist due to the relative coarse resolution of the atmospheric forcing, a potentially still too low resolution of shelf areas and the neglection of tidal forcing.

375
However, as we will see in the following sections, the more realistic circulation regime in higBUS has a significant impact on the simulated biogeochemical distributions and the transient response to anthropogenic-induced changes.

Biogeochemical pattern
The pattern of biogeochemical variables in the BUS are primarily controlled by water mass composition and upwelling structures (Hutchings et al., 1995;Williams and Follows, 2003). Thus, the pattern of net primary production (NPP) resembles the 380 spatial distribution of upwelling with high NPP rates along the coast in higBUS (Fig. 10).
In the upwelling center of Lüderitz (26-28 • S), the NPP in higBUS is significantly higher (times 1.5-2) in the near shore band (about 155 km wide, Fig. 10c) than in GR15. Off Namibia at 23 • S we also find a stronger NPP while the climatological mean upwelling rate is comparable between higBUS and GR15 at this location.
High NPP can only be achieved by high nutrient supply to the production zone. Subsequent high export rates of organic 385 material potentially lead to high remineralization induced oxygen consumption at depth. This is reflected in the low oxygen concentrations on the shelf in higBUS (Fig. 11). The offshore regions show very similar concentrations as found in the climatology from Schmidtko et al. (2017). Note, that the shelf areas of the BUS might not be well represented in the climatology from Schmidtko et al. (2017). In GR15, we find a broad band of low O 2 concentrations along the entire coast of the northern BUS (north of 30 • S), stretching out into the open ocean. Moreover, the O 2 concentrations in this 2-10 • broad band are lower 390 than the threshold value for nitrate reducing processes (Fig. 11c, black line) and even pass the critical O 2 concentration of 2 mmol m −3 for complete denitrification (Fig. 11c, lightblue line). This is a first indication that different biogeochemical processes are at play in GR15 and higBUS.
To evaluate the performance of higBUS and to illustrate the very different near shore conditions between higBUS and GR15, we focus on the well studied area off Namibia (box limits at 22 • -24 • S and 12 • E -coast, see e.g. Emeis et al. (2018)). We 395 compare model data to mean concentrations calculated from data of 4 GENUS cruises (Lahajnar (2012a(Lahajnar ( , b, 2015a) for the same region. Cruise data and model data were sampled and averaged to a 0.5 • by 0.5 • grid with vertical resolution of the model. Monthly mean model data are taken only where observational data are available.
Despite higher local NPP (+12 %) and, thus, higher export of organic material, we find higher oxygen concentrations at depth in higBUS than in GR15 (Fig. 12a). O 2 concentrations in higBUS are lower than the observations from 4 GENUS campaigns 400 (Lahajnar, 2012a(Lahajnar, , b, 2015a, but within the spatial variability. Most important, the mean oxygen concentration is above 20 mmol m −3 , which is our defined threshold value for the occurrence of nitrate reduction processes. The nitrate and nitrite profiles of higBUS fit very well to the observation, indicating a prevailing remineralization of organic material on dissolved oxygen or a fast reoxidation of products from nitrate reduction processes. Simulated rates of the first step for nitrate reduction (DN RN ) and nitrification of NH + 4 and NO − 2 fit well to observed rates from other coastal OMZs  and 405 model studies (Azhar et al., 2014). Nitrogen loss processes as denitrification (N RN 2) or anammox (Anam) mainly occur between 100-150 m depth at very low rates in higBUS. In contrast, in GR15 we see high rates of the nitrate loss processes (DN RN , N RN 2, and Anam) occurring in a very low oxygen environment (Fig. 12a) and resulting in NO − 3 concentration of only 15 % of the observed concentration at 400 m (Fig. 12c). Additionally, NO − 2 accumulates to more than the 40-fold observed value. The lack of oxygen prevents nitrification back to NO − 3 .

410
GR15 also shows the typical feature of artificial nutrient trapping, that is present in upwelling regions in many coarse resolution models indicated by accumulated phosphate concentration (Najjar et al., 1992;Aumont et al., 1999;Dietze and Loeptien, 2013). The local N:P ratio below 50 m is significantly lower in GR15 than observed or simulated in higBUS (Fig.   14). Below 200 m a reversed linear N:P relationship clearly indicates that strong N-loss processes are in full swing. Given the prescribed fixed N:P ratio for organic matter production in HAMOCC, it is obvious that the upwelling water in GR15 can 415 supply only a limited NPP due to the lack of nitrate (Fig. 10).
Nitrogen fixation, the marine process that, in general, compensates the N-loss from denitrification and anammox in the ocean,  BUS region, there is no published observational evidence for the presence of Trichodesmium  in line with 420 our simulations. Furthermore, no or very low nitrogen fixation was detected during different measuring campaigns (Staal et al., 2007;Wasmund et al., 2015). This comparison to observations gives us the confidence that the N-cycle is well represented in higBUS, based on sufficient oxygen supply to keep ocean water in the BUS region well above O 2 levels at which nitrogen loss processes start. In contrast, hydrodynamical processes in GR15, especially the water mass ventilation, are too weak to stabilize O 2 levels above our nitrate 425 reduction threshold (O 2 crit = 20 mmol m −3 ). The consequences are too low concentrations of oxygen and nitrate, higher than observed phosphate levels, and unrealistic high concentrations of nitrite.

Driving forces of O 2 supply
The physically induced contribution for different O 2 pattern of higBUS and GR15 is primarily lateral transport, especially in the zonal, onshore direction (Fig. 15). The O 2 transport in several zonal bands is towards the coast in higBUS. Around 17 • S 430 the zonal transport is away from the coast, but here, we find the largest meridional O 2 transport coming from the North induced by the poleward coastal undercurrent. All of these features are missing in GR15 due to the different horizontal circulation pattern (Fig. 5) and the depleted O 2 concentrations in a broad off shore region (Fig. 11). In addition, climatological mean annual maximum mixed layer depths (Fig. A4) in higBUS are deeper in the near shore coastal band. Therefore, areas on the shelf and in its vicinity are better ventilated than in GR15.  A closer look at the zonal section along 23 • S shows that the climatological mean transport of O 2 is clearly eastward, i.e. towards the coast in the upper 400 m in higBUS (Fig. 16a). Only in the upper most 20-30 m the wind induced circulation transports oxygen away from the coast. The eddy-induced O 2 transport east of 11.5 • E is also onto the shelf with a maximum between 50-150 m on the order of additional 10 % of the mean O 2 transport (Fig. 16b). Thus, the meso-scale structures at this latitude contribute to the ventilation of the shelf. In GR15, we find a similar pattern, but the calculated climatological mean 440 O 2 transport is much lower (Fig. 16c). This is due to overall smaller zonal velocities in GR15 (Fig. 6b), but also due to the more extended OMZ at 200-400 m (Fig. 11c) with lower absolute O 2 concentrations. Interestingly, the sign of the "eddy-induced" transport is negative (i.e. a westward O 2 transport). Thus, the temporal variability expressed by u ′ O ′ 2 rather lowers the zonal mean O 2 transport in GR15.
As already mentioned, the NPP on the shelf of higBUS is 3-4 fold higher than in GR15 (Fig. 17a). In addition, the zonal 445 transport of detritus follows the zonal circulation with a strong eastward component (Fig. 17b). Thus, there is a high turnover of organic material which is reflected in the O 2 consumption (Fig. 17e). Highest rates are found in the upper 50 m and indicate primarily the O 2 demand by grazers, but aerobic remineralization occurs throughout the upper 300 m of the water column. As seen for the station data (Fig. 12a), climatological O 2 concentrations stay well above O 2crit . Due to the high NPP, the air-sea flux of O 2 shows an outgassing almost over this entire section (Fig. 17d). An exception is only a small area close to the coast, 450 which show an invasion of atmospheric O 2 . Upwelling of nutrient rich water produces a high NPP, but this water is also cold and has, thus, a high solubility. Overall, the air-sea exchange rather plays a minor role in O 2 supply of the shelf at 23 • S in higBUS. Again, the spatial pattern in GR15 is quite different with a lower NPP (Fig. 17a) and lower detritus transport onto the shelf (Fig. 17c). Still, O 2 concentrations are depleted below O 2crit (Fig. 17f)    Because all parameter settings of the biogeochemical model are identical in both set-up, we can attribute the emerging differences in the tracer distributions to the higher horizontal grid resolution of higBUS which allows for the development of meso-scale circulation structures. Our finding on the impact of horizontal model resolution is in line with earlier studies of regional (e.g. Schmidt and Eggert (2016)) or global set-ups (Duteil et al., 2014;Buchanan and Tagliabue, 2021). Furthermore, 465 the statement of Séférian et al. (2020) that biases found in the O 2 distribution of Earth system models are independent of the model resolution can not be supported by our results. The adequate representation of meso-scale dynamics is essential not only to achieve a good O 2 distribution (Buchanan and Tagliabue, 2021), but, moreover, to avoid the occurrence of unrealistic nitrogen losses. As we will show in the next sections, also temporal tracer trends are affected.  (Schmidtko et al., 2017) (b), and GR15 (c). The rectangle in (a) indicates the area for the O2 budget given in Table 2.

O 2 trends in the BUS region over 1960 to 2009
470 The role of oxygen controlling the remineralisation pathways and, thus, the loss of bioavailable nitrogen gets addtional relevance in the context of climate change and ocean warming. Deoxygenation might enhance N-loss processes and could affect marine habitats with detrimental consequences of marine life (Schmidtko et al., 2017). A reduction of the global oceanic O 2 content of more than two per cent since the 1960s has been deduced (Schmidtko et al., 2017). The observed O 2 change is not uniform, but eastern boundary upwelling regions are most vulnerable with a potential of an intensification of upwelling 475 bringing more O 2 depleted water onto the shelf (Breitburg et al., 2018). The data compilation by Schmidtko et al. (2017) shows a dipole structure for the trend in O 2 of subsurface waters (200-400 m) in the BUS between 1960 and 2009 (Fig. 18b). In our simulation, a similar pattern is found in higBUS with slightly higher magnitudes and a small northward shift (Fig. 18a).  (Stramma et al., 2012;Ito et al., 2017). Interestingly, in higBUS we find a positive decadal trend south of 30 • S and west of 10 • E which is missing in Schmidtko et al. (2017), but is 485 also found in the data compilation of Stramma et al. (2012) (their Fig.2b). In our simulation, this positive trend pattern might be attributed the intensification of the wind stress curl over this region (Fig. A5). However, the poor data coverage, especially in the last decade (Stramma et al., 2012Schmidtko et al., 2017), creates a high uncertainty on the observed trend magnitude. Two mechanisms could provide additional oxygen at intermediate depth: 1) a decrease in O 2 consumption due to a reduced NPP and, thus, reduced particle flux to the corresponding depth or 2) an increased ventilation of this depth horizon either by vertical or lateral O 2 transport. NPP and particle fluxes into this box increase in both setups between the two time periods 500 (Table 2, Fig. A7). Consequentially, even more detritus is remineralized in the depth range of 200-400 m with an increase of 13 % and 23 % of O 2 consumption due to remineralization for higBUS and GR15, respectively. The higher O 2 demand is overcompensated by changes in lateral and vertical transport in higBUS, where the O 2 inventory of the box increases by more than 15 %. In contrast, O 2 inventory in GR15 decreases slightly due to enhanced particle flux, with a concurrent small increase in N-consumption due to remineralization (3 %) .

505
In line with the higher absolute O 2 transport shown in Fig. 15, lateral transport components in higBUS are higher than in GR15 ( Table 2). The most striking difference between higBUS and GR15 is the vertical O 2 transport at 200 m and 400 m into the box. While in higBUS an intense upwelling at 400 m carries oxygen into the box, the vertical transport at 200 m acts as a small O 2 loss (upward transport). In GR15, the box is primarily ventilated from above (across 200 m), while small mean upwelling velocities at the lower boundary play a minor role. However, the decreasing water mass age in higBUS over the 510 decades (Table 2) further indicates a supply of water being recently in contact with the atmosphere.
In both set-ups, the net O 2 transport into the box increases over the decades. Because the chosen box covers a region with diverging horizontal currents, the net zonal transport out of the box increases, but this is compensated by mean net meridional and vertical transports into the box. The transport changes are a result of changes in the forcing which lead to a change in the barotropic streamfunction (Fig. A6). The change in the barotropic streamfunction is also present in GR15 (Fig. A6d), but with 515 a different pattern and largest changes rather at the western side of the investigated area.
The change in the dye tracers (Fig. A8) for the region of the box indicates that the additional oxygen comes primarily from the Indian ocean (SIO) in higBUS. The intensification of the upwelling on the shelf between 10 • S to 30 • S and the off-shore westward transport north at 20 • S are clearly reflected in the change of the SIO dye tracer in higBUS, while this change is absent in GR15. The increased fraction of water from the SIO is linked to a stronger leakage of the Agulhas Current into the South 520 Atlantic at the end of the simulation period, which has also been identified by Biastoch et al. (2009). These authors attributed the changed leakage strength to a relocation of southern hemisphere westerlies between 1970 and 2009. The O 2 transport for the later period (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) and the decadal change of the O 2 transport are in line with a stronger Agulhas leakage (Fig. A9).
Both model resolutions show an enhanced O 2 transport with the Benguela current. However, we find that only in higBUS this higBUS GR15 1960GR15 -1979GR15 1990GR15 -2009GR15 1960GR15 -1979GR15 1990GR15 -2009 O2 inventory (10 10 kmol)  Table 2. Mean O2 inventory of a box region (15-20 • S and 10 • W-5 • E) over the 200-400 m depth range, zonal, meridional, and vertical mean O2 transport across the box boundaries as well as the total net transport into the box, biological induced O2 fluxes, and mean water mass properties. Positive numbers for transport indicate eastward (zonal), northward (meridional) or upward (vertical) fluxes. Bold numbers of transport are given to identify fluxes into the box more easily (e.g the negative vertical transport across 200 m in GR15 is an O2 gain for the box). Note, that the particle flux is given at 100 m.
additional O 2 reaches the box region and compensates the O 2 demand from an increased NPP.

525
Within the model framework, we conclude that changes in O 2 in the northern BUS region, which are found in the observation over the last decades, are related to changes in the water mass composition with a higher contribution of SIO water. Changes in the wind stress forcing modify the upwelling pattern and lead to higher NPP along the coast. This trend pattern is only captured by the high resolution model version. Our results underline the importance of a good representation of the physical conditions 530 for the biogeochemistry. The coarse resolution setup is not suitable to depict the mean distribution of biogeochemical tracers and to resolve their trends in this eastern boundary upwelling region.

Surface pCO 2 and decadal trends of the carbon inventory in the BUS
Recent publications raised the question, if the BUS region is rather a source or sink for atmospheric CO 2 Laruelle et al., 2018). The partial pressure difference of CO 2 between the surface ocean and the atmosphere and, thus, the 535 air-sea carbon fluxes are determined by the interplay between physical and biogeochemical conditions. With our setup we are able to investigate changes in the regional carbon cycle because we achieved a nearly steady state of our globally consistent model before we applied the transient ERA20C forcing. All horizontal and vertical tracer gradients in the global ocean are a result of the implemented physical and biogeochemical processes.
In higBUS, the entire region north of 31 • S is a source for atmospheric CO 2 for the time period 2000-2009 (Fig. 19a).

540
Upwelling of carbon rich water leads to outgassing despite the simulated high biological production rates which rather tend to reduce surface pCO 2 . The comparison to GR15 underlines this impact of biological production. While physical parameters such as coastal upwelling ( Fig. 8 and Fig. 10c, black line) and sea surface temperature of the adjacent open ocean (Fig. 9) are very similar in higBUS and GR15, the integrated NPP is higher in higBUS (Fig. 10). Thus, biological activity sequesters dissolved inorganic carbon and exports it to greater depth and thereby reduces local outgassing of upwelled carbon rich water.

545
The difference in NPP is related to the lack of bioavailable nitrogen in GR15, as has already been shown (Fig. 14), and which is a consequence of low oxygen conditions (see also Fig. A10 of the N/P ratio at 100 m). In higBUS, high ∆pCO 2 values (> 20 ppmV, ∆pCO 2 = pCO 2 in seawater minus pCO 2 in air) are found closer to the coast in a 400-500 km wide band. In GR15, higher positive ∆pCO 2 values also stretch out far to the west which is not confirmed by measurements of previous pCO 2 compilations for the global open ocean (Takahashi et al., 2009;Bakker et al., 2016). South of 32 • S both setups show a very 550 similar small negative ∆pCO 2 . The annual mean carbon flux is shown in the Appendix for completeness (Fig. A12).  In this study, we address the temporal evolution of the biogeochemistry in the Benguela upwelling system (BUS) over the 20th century. Therefore, we set up a global ocean biogeochemical model with a stretched grid configuration which allows also for resolving meso-scale circulation structures in the South Atlantic (higBUS). The advantages of this set-up are manifold: a) it captures the Agulhas leakage as well as the impact of eddies which shed off the Agulhas current (Biastoch et al., 2009;590 Jungclaus et al., 2013) and carry oxygen rich water from the Indian Ocean into the BUS; b) it represents the complex frontal system and the poleward undercurrent of the BUS (Monteiro et al., 2011;Gutknecht et al., 2013); and c) its global approach avoids artefacts due to prescribed oceanic boundary conditions and allows for long and consistent model simulations, especially important in view of the marine carbon cycle. We set an additional focus on the development of the oxygen minimum zone (OMZ) and nitrogen loss processes within the BUS. Therefore, we include a prognostic treatment of ammonium and nitrite and 595 associated conversion processes like nitrification and anammox. The impact of meso-scale structures on the biogeochemistry is identified by running the same model in a coarser set-up (1.5x1.5 • , GR15) where the impact of eddies is parameterized with the common scheme from Gent et al. (1995).
The evaluation of higBUS with climatological data, cruise data, and satellite derived products shows a high level of agreement. In particular, the prognostic tracers NO − 2 and NH + 4 provide an additional assessment option for the role of anaerobic The impact of meso-scale structures, i.e. eddies and filaments, enhances vertical mixing and leads to a ventilation of the upper ocean. In higBUS, we find low O 2 (< 20 mmol/m 3 ) only on the shelf, which is in agreement with observations. Furthermore, the large scale circulation pattern and water mass properties in higBUS differ form those in GR15. The consequences of a different mean circulation are evident in the multi-decadal trends of the oxygen. Only in higBUS, we find an 610 increase in mid-depth O 2 being also reported in observational trend estimates (Schmidtko et al., 2017;Stramma et al., 2012;Ito et al., 2017). The results of higBUS are in agreement with model studies of age and pathways of the water masses in BUS of Tim et al. (2018) who applied water parcel tracking methods over the period 1960-2009. They speculated on biogeochemical changes such as increased oxygen supply of the BUS induced by the intensification of the Agulhas leakage which is supported by our findings from higBUS.

615
All these resolution dependent differences eventually affect the carbon cycle and the water column carbon inventory change over the 20th century. In principle, we find for both model resolutions that the northern BUS is a carbon source for the atmosphere, while the southern BUS is rather a sink. However, the source region in GR15 is much more extended into the South Atlantic, which disagrees, at least for the region south of 22 • S, with the observational estimates (Takahashi et al., Despite a very similar global uptake of anthropogenic carbon in both set-ups, the regional C-inventory increase over the entire simulation is smaller in higBUS. For the last 5 decades ) higBUS projects a negative trend close to the coast.
The different trends of the carbon chemistry in higBUS and GR15 also show up in the local pH trend. In higBUS, pH changes on the shelf north of 32 • S are even slightly positive. Thus, according to our model, there is no impact of ocean acidification over the last decades in the BUS. Unfortunately, there is no long term data set on pH for this region to evaluate our finding.

625
Despite the better performance of the set-up with meso-scale structures, we have to consider that the higBUS experiment has several shortcomings: (1) The spatial resolution of the chosen forcing from reanalysis data (ERA20C) is rather coarse (1 degree) which might have the tendency to underestimate coastal upwelling as discussed by Milinski et al. (2016). Especially, a wind product with higher resolution, i.e. QuikSCAT and ASCAT with 0.25 • resolution, might enhance the upwelling strength (Junker et al., 2015). Unfortunately, these higher resolved data products cover only one decade and are, thus, not suitable for our 630 purpose.
(2) The biogeochemical model includes a comprehensive representation of the nitrogen cycle, but it does not include remineralization processes based on the activity of chemolithoautothrophic sulphur bacteria (see e.g. Schmidt and Eggert (2016)). They are active in anoxic sediments layers and could supply ammonium to the bottom layer from DNRA (dissimilatory nitrate reduction to ammonium).
(3) Furthermore, our setup does not include migrating zooplankton. Zooplankton degrades organic material and consumes O 2 . Migration could therefore have an impact on the oxygen distribution in the water column.

635
However, we are not convinced that the behavior of migrating zooplankton in the vicinity of OMZs is well enough understood, yet (Escribano et al., 2009).
In summary, our results emphasize that it is crucial to resolve meso-scale circulation structures to achieve a reasonable oxygen distribution and, concurrently, a good representation of the nutrient cycle in the BUS. The coarse resolution version, GR15, suffers from too little mixing and too low exchange between the open ocean and the shelf areas. The resulting extended 640 OMZ affects not only the nitrate reduction processes, but also the local primary production and, thus, the biogeochemistry of the BUS. By simulating a comprehensive nitrogen cycle these weaknesses of the coarse resolution set-up clearly shows up in the pattern of NO − 2 . Thus, in view of a potentially increasing number of NO − 2 measurements, a comprehensive nitrogen cycle could be a promising validation tool for biogeochemical models. In addition, we found that the impact of meso-scale structures also changes the transient behavior of the local O 2 distribution, the regional CO 2 uptake, and the temporal change of 645 sea water pH. Therefore, we are convinced that, at least for some regions, the projected future pathway of changes in the ocean biogeochemistry cannot be reliably derived from coarse resolution models. given, we adjusted the number, otherwise we take the value from the given reference. (a) Lima and Doney (2004) We use data from four GENUS cruises:

A3 Climatological pattern of water mass composition
Fig . A1 shows the water mass composition at 30 • S for higBUS and GR15 for water from 6 regions. Water from SIO has the largest contribution on the shelf and off shore. Only the dye tracer (sBUS), which origins in the area that hosts this section, has higher values in the upper water column. Water coming from SATL contributes to 15-20 % to shelf water (<400 m depth). In 660 higBUS, there is also some fraction ( < 20 %) of water from SACW reaching 30 • S on the shelf.  A2 shows the water mass composition at 23 • S for higBUS and GR15 for water from 6 regions. As shown in Fig. 6, water from the SIO contributes the largest fraction west of 12 • E below 150 m. This is a result of the transport by the Benguela current which is clearly identifiable in higBUS, but has a more diffusive structure in GR15. About 20 % of the water mass origins from SATL, being mixed in the Benguela current west of the Cape Basin (Fig. A1) in both resolutions. In higBUS, the 665 shelf (< 400 m) is primarily occupied by water from SACW with additional small contributions of water from the equatorial Atlantic (EATL, 5 %) and nBUS (>10 %) at this northern BUS section . The dye tracer, that origins in the BUS part which also hosts this section (sBUS), shows, of course, the highest fraction in the upper 100 m of the water column. However, the more efficient upwelling and the concurrent stronger westward transport in higBUS is reflected in the pattern of this dye tracer, which is confinded to the upper 80 m of water column on the shelf. In contrast, in GR15 we find a more diffusive picture with 670 a higher fraction of this dye tracers found in the bottom water on the shelf.

A4 Presence of eddy structures in higBUS
The impact of the higher horizontal resolution on the flow field becomes visible by the sea surface height (SSH) anomaly, here shown for the southern BUS (Fig. A3). In this daily mean snapshot for April, 10th 2009 higBUS shows meso-scale eddy structures of 100 km diameter which are travelling westward into the Atlantic. The clear displacement of the anticyclonic 675 eddies is highlighted by the 3 cm contour line of the SSH-anomaly also shown for 10 and 20 days after April, 10th. In GR15, we see a much broader positive SSH anomaly which is rather stationary and fully flattens within 20 days (see only one contour line in Fig. A3b). As seen for surface (Fig. 7a), subsurface (Fig. 5a) and upwelling velocities (Fig. 8) higBUS also shows a SSH anomaly pattern with higher variability indicating a vivid mixing of water masses, especially south of 30 • S. Note non-linear color scale. A7 Temporal change of the NPP and the water mass distribution Figure A7. Temporal change in NPP (gC m −2 d −1 ) between (1990-2009) and (1960-1979) for higBUS (c) and GR15 (d).
As shown for the zonal section of 23 • S in Fig. A2, water with the origin from SATL contributes 20-25 % to the water mass composition of the BUS in the northern part (north of 30 • S) while the southern BUS gets less than 10 % of this water at a mean depth of 100-200 m in higBUS (Fig. A8). Water from SACW arrives with the shelf-bound poleward current into the 685 BUS. Therefore higher fractions are found only on the shelf (Fig. A2) in higBUS. The temporal change of both dye tracers, SATL and SACW, is different from the change we found in O 2 which underlines the finding that it is primarily the change of the Agulhas leakage that is responsible for the O 2 anomaly north of 20 • S. The large scale patterns of both tracers in GR15 have some similarity to higBUS with the already discussed deviations, but the temporal change is different. Particularly, water coming from the SATL has a slightly higher impact on coastal regions south of 30 • S in the later period 1990-2009. This might 690 explain the increased NPP south of 30 • S (+20 %, Fig. A7) which causes also a pronounced O 2 change not being observed at that magnitude.    [1960][1961][1962][1963][1964][1965][1966][1967][1968][1969][1970][1971][1972][1973][1974][1975][1976][1977][1978][1979] for the same dye tracers: SATL (g,j), SIO (h,k), and SACW (i,l) for higBUS (g,h,i) and GR15 (j,k,l). Contour line is given at 0. Figure A9. As Fig. 15, but absolute mean lateral O2 transport averaged over 1990-2009 (a,b) and the difference of the absolute mean lateral O2 transport between (1990-2009) and (1960-1979) (c,d) for higBUS (a,c) and GR15 (b,d).

A8 N/P ratio and climatological remineralization of detritus
The comparison of the N/P ratio at 100 m clearly shows the depleted coastal area in GR15 (Fig. A10). There are some discrepancies in the pattern between higBUS and the WOA data. However, the ratio is very sensitive to quality of phosphate 695 measurements. Another N/P estimates from data miss e.g. the very low ratio west of 12 • E (Galbraith and Martiny, 2015).