A Novel Modeling Approach to Quantify the Influence of Nitrogen Inputs on the Oxygen Dynamics of the North Sea

1. Introduction

The availability of dissolved oxygen (O₂) in marine waters is crucial for the metabolism of many marine species (e.g., Rosenberg et al., 1991; Díaz and Rosenberg, 1995). Hypoxic conditions, i.e., dissolved O₂ concentration < 2 mg O₂ L⁻¹ (e.g., Tyson and Pearson, 1991; Rabalais et al., 2002), constitute a major stressor to marine species and can have a severe impact on an ecosystem's biodiversity (Gray et al., 2002; Vaquer-Sunyer and Duarte, 2008; Topcu et al., 2009; Friedrich et al., 2014).

In the North Sea — a northwestern European shelf sea — hypoxic conditions occurred in the summers of the 1980s (Rachor and Albrecht, 1983; von Westernhagen and Dethlefsen, 1983; Brockmann and Eberlein, 1986; Peeters et al., 1995), causing death of benthic and demersal species (von Westernhagen and Dethlefsen, 1983). Different studies showed that these occurrences of low O₂ conditions evolved during seasonal stratification (von Westernhagen and Dethlefsen, 1983; Peeters et al., 1995), driven by eutrophication due to excess riverine nutrient loads (Brockmann and Eberlein, 1986; Brockmann et al., 1988; Peeters et al., 1995). In consequence, the 2nd International Conference on the Protection of the North Sea declared to “aim to achieve a substantial reduction (of the order of 50%) in inputs of phosphorus (P) and nitrogen (N) to these areas between 1985 and 1995” (ICNS-2, 1988), in order to mitigate the effects of eutrophication in the North Sea.

Despite significant reductions in riverine nutrients, especially for P but not for N (Artioli et al., 2008; Claussen et al., 2009; OSPAR, 2013), seasonal O₂ deficiency (i.e., dissolved O₂ concentration < 6 mg O₂ L⁻¹; OSPAR, 2003) still occurs regularly in the southeastern North Sea (Topcu and Brockmann, 2015). Ongoing O₂ monitoring even suggests recurrent (close-to-)hypoxic conditions (Weigelt-Krenz, pers. comm.), though, scientific publications on such events are sparse (Weston et al., 2008). Thus, O₂ deficiency remains a persistent problem in the North Sea and the investigation of its causes is essential for the ecological management.

Große et al. (2016) applied a physical-biogeochemical model to the North Sea analyzing the physico-biochemical interactions that make the southeastern North Sea most susceptible to O₂ deficiency. They found that the occurrence of distinct seasonal thermal stratification, resulting in a thin sub-thermocline layer in this region of intermediate depth, and high local near-surface net primary production (NPP) constitute the key factors. Consequently, near-bottom O₂ concentrations can become especially low during years of very high NPP. High NPP enhances (pelagic and benthic) bacterial re-mineralization of organic matter, consistently accounting for more than 80% of near-bottom O₂ consumption in the region of lowest bottom O₂ concentrations (about 55.5–56.0 °N, 6–7 °E).

In contrast to near-shore eutrophication problem areas affected by high river loads entering the coastal areas, the off-shore areas affected by O₂ deficiency are exposed to a wider variety of nutrient sources. Thus, the key question is which nutrient inputs drive NPP and O₂ consumption finally, resulting in O₂ deficiency in these North Sea regions?

Past studies suggested that nutrient inputs from adjacent North Sea regions may have contributed to the evolution of hypoxia in the 1980s (Brockmann and Eberlein, 1986; Brockmann et al., 1988). More recent publications also pointed on nutrient inputs from adjacent seas, especially the Northeast Atlantic (Artioli et al., 2008; Gröger et al., 2013). In addition, the model study by Troost et al. (2013) highlighted the influence of atmospheric N deposition on NPP in the southeastern North Sea and on near-bottom O₂ concentrations. Quantifying the contributions of these different nutrient sources in the regions susceptible to O₂ deficiency is vital for the development of potential management strategies for these regions. This demand is in line with the “integrated target-oriented and source-oriented” management approach by OSPAR (1999).

A nutrient tagging method applicable to biogeochemical models — often referred to as “trans-boundary nutrient transports” (TBNT; e.g., Blauw et al., 2006; OSPAR, 2010) — facilitates the tracing of elements (e.g., N) from individual sources throughout the entire biochemical cycle and the entire ecosystem (Ménesguen et al., 2006). By this, the TBNT method enables the quantification of the influence of these individual nutrient sources on biochemical quantities (e.g., total nitrogen; TN) and processes (e.g., NPP) in different regions of an ecosystem.

This method has been applied to various European marine ecosystems, including the Baltic Sea (e.g., Neumann, 2007; Radtke et al., 2012) and the North Sea (e.g., Blauw et al., 2006; Lacroix et al., 2007; Painting et al., 2013; Troost et al., 2013). Troost et al. (2013) demonstrated the importance of atmospheric N deposition for the off-shore regions of the southern North Sea, where it contributed 10–20% to overall TN and even 20–30% to NPP in 2002. Other studies (OSPAR, 2010; Painting et al., 2013) showed that N inputs from the large Dutch and German rivers are the major contributors to TN in the coastal southern and southeastern North Sea, respectively, while those from the rivers at the British North Sea coast extend far east into the off-shore southern North Sea. This, together with the findings by Große et al. (2016), suggests that these N sources have a strong influence on the O₂ dynamics in the southern and southeastern North Sea.

Up to now, the TBNT method was limited to the tracing of elements brought into the system by the selected sources, not allowing for the consideration of O₂. Thus, quantitative analyses of the influence of different nutrient sources on the O₂ dynamics required additional model scenarios within which individual nutrient sources were removed (e.g., Troost et al., 2013) or nutrient loads were reduced (e.g., Lenhart et al., 2010). This implies an individual scenario was required for each individual nutrient source in order to assess its impact on the O₂ conditions which is time-consuming and causes high computational costs.

Therefore, this study presents a new approach for expanding the TBNT method by a direct link of the O₂ dynamics to the nutrient inputs from individual sources. We focus on the relative importance of riverine vs. non-riverine sources of N, and their influence on the O₂ dynamics in the regions most susceptible to O₂ deficiency.

We first describe the applied HAMSOM-ECOHAM model and provide the relevant information on the theory and the implementation of TBNT. Große et al. (2016) showed that the model system captures the main features of observed O₂ concentrations in the North Sea. Therefore, only the simulated volume transports are validated and spatial distributions of TN from selected sources are presented and compared to earlier TBNT studies. Finally, the results of the expanded TBNT method are presented with special focus on two regions regularly affected by O₂ deficiency — the southeastern North Sea (Topcu and Brockmann, 2015) and the Oyster Grounds (e.g., Greenwood et al., 2010).

2. Materials and Methods

2.1. The HAMSOM-ECOHAM Model

We use a three-dimensional (3D) model consisting of the physical model HAMSOM (HAMburg Shelf Ocean Model; Backhaus, 1985; Pohlmann, 1991, 1996) and the biogeochemical model ECOHAM (ECOsystem model, HAMburg; Pätsch and Kühn, 2008; Kühn et al., 2010; Lorkowski et al., 2012; Große et al., 2016). HAMSOM is a baroclinic, primitive equation model using the hydrostatic and Boussinesq approximations (Pohlmann, 1991). The current velocities are calculated using the component-upstream scheme. The horizontal is discretized on a staggered Arakawa C-grid (Arakawa and Lamb, 1977) and z-coordinates are applied to the vertical. A detailed description of HAMSOM is provided by Pohlmann (1991, 1996, 2006). HAMSOM applications can be found, e.g., in Backhaus and Hainbucher (1987), Chen et al. (2013), Pohlmann (1996), and Mathis and Pohlmann (2014). The 3D fields of temperature (T), salinity (S), advective flow and vertical turbulent mixing coefficients calculated by HAMSOM are used as forcing for ECOHAM.

The ECOHAM model represents all parameter groups of a so-called NPZD-type model (nutrients-phytoplankton-zooplankton-detritus) to represent the lower trophic level dynamics. It describes two phytoplankton groups (diatoms, flagellates) and two zooplankton groups (micro- and mesozooplankton). In addition, the “microbial loop” (Azam et al., 1983) is represented. ECOHAM simulates the cycles of carbon (C), N, P, silicon (Si) and O₂ (Pätsch and Kühn, 2008; Lorkowski et al., 2012). The present version of ECOHAM is identical to that used by Große et al. (2016). Lorkowski et al. (2012) provide a full description of the ECOHAM model equations and parameter settings. Here, we will only summarize the relevant features of the O₂ cycle.

Figure 1 illustrates the O₂ cycle implemented to ECOHAM (Müller, 2008; Pätsch and Kühn, 2008), including the links to the cycles of C, N and P. The arrows indicate the loss or gain of a substance, e.g., $\text{N}{\text{H}}_4^{+}$ is released during benthic re-mineralization and utilized during pelagic nitrification. The air-sea flux of O₂ is parametrized according to Wanninkhof (1992). In the interior water column, NPP according to Liebig (1841) and Steele (1962) constitutes a source process for O₂ and links O₂ to the cycles of N, P, and C. All other biochemical processes, namely respiration by zooplankton (Sterner and Elser, 2002) and aerobic bacteria, and nitrification, represent O₂-consuming processes. While the former processes constitute links of O₂ to the cycles of C, N, and P, nitrification is only linked to N. Pelagic denitrification is implemented to ECOHAM, but is negligible as it only takes place under anaerobic conditions, not occurring in the North Sea. Pelagic anaerobic $\text{N}{\text{H}}_4^{+}$ oxidation (anammox) is not implemented for the same reason (Große et al., 2016). NPP is T-dependent and pelagic nitrification is light-limited in the applied model setup.

The sediment is described by a simple zero-dimensional module (Pätsch and Kühn, 2008). Benthic re-mineralization follows a first-order approach inhibiting year-to-year accumulation of matter (Große et al., 2016). Benthic denitrification is linked to benthic O₂ consumption according to Seitzinger and Giblin (1996), reducing the O₂ concentration in the deepest pelagic layer (Fennel et al., 2006). Explicit benthic nitrification and benthic anammox are not implemented (Pätsch and Kühn, 2008).

In relation to the re-mineralization processes it has to be mentioned that the C:N:P ratios of particulate (POM) and dissolved organic matter (DOM) can evolve freely with time. POM is represented by two detritus species (small/large). Model DOM and POM only represent the reactive fractions (Kühn et al., 2010). Carbon oxidation occurs with a molar ratio of C/O₂ = 1 (Neumann, 2000; Paulmier et al., 2009). Therefore, the O₂ consumption due to respiratory processes can be described by the corresponding net C release, except for nitrification for which the molar ratio between $\text{N}{\text{H}}_4^{+}$ uptake and O₂ consumption is 2. The model equations and parameters of the biochemical processes affecting O₂ are given in Appendix A (Supplementary Material).

2.2. The TBNT Method

The basic concept of TBNT is that any property (e.g., input source) of an element (e.g., N) can be traced throughout the entire elemental cycle represented by a biogeochemical model. This is done by introducing additional model state variables and processes for each property. This concept dates back to the work by Ménesguen and Hoch (1997), who introduced the age as a property, and later on has been applied to various models and ecosystems (e.g., Wijsman et al., 2004; Blauw et al., 2006; Ménesguen et al., 2006; Lacroix et al., 2007; Neumann, 2007; Timmermann et al., 2010; Radtke et al., 2012; Painting et al., 2013; Troost et al., 2013; Perrot et al., 2014; Radtke and Maar, 2016).

According to Ménesguen and Hoch (1997), any property p (e.g., origin/source etc.) can be attached to any state variable X, and can be traced within a biogeochemical model by solving an additional differential equation for the product X · p instead of the bulk state variable X. The product X · p then represents that subset of the state variable X with the defined property p. Technically, the product X · p is introduced into the model as a new fraction state variable X_p = X · p and the related processes are introduced according to the differential equation describing the changes in the corresponding bulk state variable X. Each combination of state variables X and properties p_i (e.g., different input sources), $X_p^i=X·p_i$, requires an additional differential equation to be solved by the model. For easier readability, $X_p^i$ is hereafter denoted as Xⁱ.

Compared to existing TBNT studies, we apply a slightly different approach which allows for the diagnostic calculation of the temporal evolution and spatial dispersal of X_i. In the following, the theory for this element tracing approach is briefly described. In this context, the terms “bulk” and “fraction” are used to refer to the overall state variables and processes, and their labeled counterparts, respectively. Thereafter, the approach for linking the biochemical processes affecting O₂ to the cycle of a traced element is presented.

2.2.1. Tracing Elements from Individual Sources

The convection-diffusion equation describes the change in concentration of a selected bulk state variable with time due to different physical and biochemical processes. This equation represents the differential equation solved by physical-biogeochemical models. A detailed derivation of the convection-diffusion equation for a fraction variable starting from that for a bulk state variable is provided by Radtke (2012). The source-specific convection-diffusion equation for a constant volume V as used in this study reads as follows:

Here, $C_X^i$ and C_X represent the concentrations of the fraction of state variable X originating from the i-th input source and that of the corresponding bulk state variable, respectively. The diffusive transport is calculated according to Fick's first law, with the second order diffusion tensor (or diffusivity) $\overline{\overline{D}}$. In the convective/advective transport term, $\overrightarrow{\nu }$ represents the three-dimensional velocity vector. R_CX represents the change in concentration of X due to the sources and sinks (i.e., biochemical processes, input from external sources). The index X_in in the fraction of this term indicates that the relative contribution of the state variable, that is consumed by a biochemical process, is used.

For n individually labeled input sources (i.e., i = 1, 2, …, n – 1, n), the concentration of each bulk state variable X at a specific location and point in time equals the sum of its contributing fraction state variables:

According to Equation (1) all physical (diffusion and transport) and biochemical processes (sources/sinks) affecting $C_X^i$ are calculated as the product of the process based on C_X and the relative contribution of the fraction variable, $C_X^i/C_X$. Thus, the major assumption of Equation (1) is that there is no preference for any fraction variable Xⁱ by any (physical or biochemical) process, which relates to the fact that all fraction variables Xⁱ are chemically and physically identical.

Compared to previous TBNT studies (e.g., Ménesguen et al., 2006; Radtke et al., 2012) this implies a different treatment of the diffusive transport. In our study its calculation bases on the spatial gradients of the bulk variable concentration. In earlier studies the concentration gradients of the fraction variable were used. In fact, it is difficult to state which approach is more correct than the other as the diffusive transport according to Fick's first law is a net transport, i.e., the actual gross transports are not known. However, it has to be noted that the diffusive transports of the individual fraction variables Xⁱ according to Ménesguen et al. (2006) and Radtke et al. (2012) can be opposed, while those following Equation (1) are unidirectional. Despite this different treatment the resulting differences in the fraction variable distributions are expected to be generally small and noticeable differences may only occur in the case of large differences between the spatial gradients of the fraction and bulk variables. The diffusive transport of the bulk variable will be the same in both cases.

The advantage of the calculation based on Equation (1) is that the element tracing can be conducted by means of a post-processing software utilizing the standard output of the applied model. It only requires the provision of the 3D output fields of all processes and state variables involved in the selected elemental cycle.

2.2.2. Linking Oxygen-Affecting Processes to Labeled Nutrients

In order to allow for the quantification of the influence of (organic and inorganic) nutrients from individual sources on the O₂ dynamics the establishment of a link between the two is required. This implies that only the local biochemical processes affecting O₂ can be linked to labeled nutrients, while physical processes (i.e., lateral effects) that only depend on O₂ itself cannot be considered explicitly.

According to Equation (1) biochemical processes on fraction variables (e.g., uptake of $\text{N}{\text{O}}_3^{-}$ containing N from the Elbe by diatoms) can be calculated using the bulk process and the relative contribution of the fraction variable to the bulk variable:

Similarly, Equation (3) can be used for linking the O₂-affecting processes to labeled nutrients. In case of ECOHAM, the O₂ dynamics are linked to the cycles of C, N, and P by different processes (see Figure 1). In case of O₂-affecting processes whose formulation is based on a state variable containing the labeled element, this state variable can be used for the link to the O₂ cycle. For instance, O₂ consumption due to nitrification is based on $\text{N}{\text{H}}_4^{+}$ containing the labeled element in case of N labeling.

For processes whose formulation is based on a different element than the labeled one (e.g., NPP with ‘currency’ C instead of N), a different approach is required. In ECOHAM, the formulation of the change in O₂ due to these processes is based on quantities that contain the currency and the labeled element, e.g., diatoms containing C, N, and P. Assuming that the concentration of the currency of the process is proportional to the currency of the labeled element these model quantities can be described as “mediators,” X_M. These mediators can be used for establishing the link between the labeled element and the O₂ cycle. For O₂-affecting processes linked to these labeled mediator variables Equation (3) then reads as:

with the O₂-affecting, biochemical bulk process R_CO2, and $C_{X_M}^i$ and C_XM representing the concentrations of the mediator fraction and bulk variable, respectively.

Using Equation (4), the following mediator variables are used for the O₂-affecting processes shown in Figure 1: diatoms and flagellates for NPP, micro- and mesozooplankton for zooplankton respiration. For pelagic re-mineralization (i.e., bacterial respiration), the bacteria are used while for benthic re-mineralization the amount of detrital matter in the sediment is used. In case of N labeling in ECOHAM, the mediator state variables are: diatom-N and flagellate-N for NPP, micro- and mesozooplankton-N for zooplankton respiration, bacteria-N for pelagic re-mineralization and sediment-N for benthic re-mineralization. In case of P labeling, which is not part of this study, the corresponding model state variables can be used as mediator.

The use of the mediator implies that the concentration of the C-based mediator variable (e.g., diatom-C) is proportional to the corresponding concentration of the N-based mediator variable (e.g., diatom-N), according to the applied stoichiometry. Conceptually, their use can be interpreted such that the biomasses (containing C, N, and P) of the different organisms (phyto-/zooplankton, bacteria) and that of benthic detrital matter are controlled by N from the individual input sources.

Our study focuses on the influence of N on near-bottom O₂ in the off-shore regions of the North Sea, where NPP is negligible due to light limitation (Große et al., 2016). Therefore, we use gross O₂ consumption (GOC), i.e., the sum of all biochemical O₂-consuming processes as the key quantity in our analysis. In ECOHAM, these processes are: pelagic and benthic re-mineralization, zooplankton respiration and nitrification. The here used GOC is comparable to the often used “apparent oxygen utilization (AOU)”.

2.3. Model Setup and TBNT Calculation

2.3.1. Model Setup

The model setup used for this study is identical to that described in Große et al. (2016). Both, HAMSOM and ECOHAM, are applied to the same model domain and grid, comprising the entire North Sea, large parts of the northwest European continental shelf and parts of the adjacent Northeast Atlantic Ocean (see Figure 2). The domain extends from 15.250 °W to 14.083 °E and from 47.583 °N to 63.983 °N. The horizontal resolution is 1/5 ° (1/3 °) with 82 (88) grid points in latitudinal (longitudinal) direction. The vertical dimension is resolved by 31 z-layers with a surface layer of 10 m average thickness. Below the surface layer, layer thicknesses increase from 5 m to 1,000 m in the deepest layers.

HAMSOM uses monthly, climatological distributions of T and S based on the World Ocean Atlas 2001 (Conkright et al., 2002) for initialization and open boundaries. The open boundary treatment for advective flow as well as T and S is described in Chen et al. (2013). The meteorological forcing is derived from NCEP/NCAR reanalysis data (Kalnay et al., 1996; Kistler et al., 2001) and provides 6-hourly information on air temperature, cloud coverage, relative humidity and wind velocity components.

For ECOHAM, short-wave radiation is calculated from astronomical radiation and cloud coverage applying a correction factor of 0.9 to account for underestimated cloud coverage in the NCEP data in temperate zones. For the biogeochemical state variables a climatology of depth-dependent monthly averages was prescribed at the boundaries and annually varying data was provided for dissolved inorganic carbon (DIC; Lorkowski et al., 2012). A daily climatology of suspended particulate matter was derived from results of Pohlmann and Puls (1994) to include its effect on the light climate.

Daily freshwater run-off and nutrient load data for 254 rivers were provided by Sonja van Leeuwen (pers. comm.) and constitute an updated dataset of that used by Lenhart et al. (2010) covering the period 1977–2014. River inputs of particulate organic N (PON) were calculated as TN minus $\text{N}{\text{O}}_3^{-}$ and $\text{N}{\text{H}}_4^{+}$, which are all available in the original data. Dissolved organic N is considered to be zero. The same approach was applied for organic P loads. Particulate organic C loads were derived from the PON loads using a C:N ratio of 106:16 (Redfield, 1934). A detailed description of the nutrient load calculation for the German and Dutch rivers is given in Pätsch et al. (2016).

The atmospheric N deposition was derived following Große et al. (2016), using annual data from the EMEP (Cooperative program for monitoring and evaluation of the long-range transmissions of air pollutants in Europe) model covering the period 1995–2012 and long-term trends (Schöpp et al., 2003). Atmospheric deposition is implemented as inputs of $\text{N}{\text{O}}_3^{-}$ and $\text{N}{\text{H}}_4^{+}$.

Both model compartments ran over the period 1977–2014. A time step of 10 minutes was applied for the HAMSOM simulation and daily fields of T, S, advective flow and the vertical turbulent diffusion coefficient were stored as output. This daily output was used as physical forcing for ECOHAM, which ran off-line applying a time step of 30 minutes. Daily ECOHAM output was stored for the entire N cycle and the biochemical processes affecting O₂.

2.3.2. TBNT Setup

For the TBNT calculation, a sub-domain within the model domain was defined, comprising the North Sea, the English Channel and the Kattegat/Belts. The outer boundaries of this sub-domain (hereafter “North Sea domain”) are: North Atlantic (NA) in the North, English Channel (EC) in the Southwest and Baltic Sea (BS) in the East (see Figure 2). An individual input source was defined for each of these outer boundaries and matter entering the North Sea across these boundaries is labeled accordingly. One atmospheric source for the entire North Sea domain was defined to include atmospheric N deposition. The individual rivers inside the North Sea domain were grouped to 8 river groups according to OSPAR (ICG-EMO, 2009; OSPAR, 2010; Los et al., 2014). One additional river group was defined collecting all rivers not attributed to one of the OSPAR groups (hereafter “other rivers/others”). The input locations of the rivers of the individual groups are indicated in Figure 2 (see different markers). Table 1 provides a complete list of all source groups.

The TBNT analysis was conducted for the period 2000–2014 using the daily ECOHAM output. A sensitivity study with an ECOHAM output time step of two hours showed that the effect of the output time step on the TBNT results is negligible (not shown). All labeled variables (except for the other rivers) are initialized with a relative contribution of 0% to the corresponding model state variables. The other rivers are initialized with a relative contribution of 100% within the entire North Sea domain. This implies that the matter initially present inside the domain is fully attributed to the other rivers, while the matter originating from the specifically defined sources has to enter the system successively.

Starting from the initial distribution, the year 1999 was re-run seven times with the final distribution of the relative contributions of each labeled state variable used as initialization for the subsequent one. This implies a minor shift in the actual masses of individual labeled state variables between the end and the beginning of two subsequent spin-up years. After the seventh re-run a quasi-steady state was reached. This gives confidence (1) that realistic distributions of the fraction variables are reached and (2) that the matter attributed to the other rivers actually relates to these rivers and not to the initial mass. The advantage of this procedure is that the TBNT analysis can be applied to a comparably short model simulation requiring only one additional year prior to the analysis period. A comparison of a spin-up over a sequence of seven years showed negligible differences (not shown) justifying this procedure.

The TBNT calculation utilizes a base time step of one day (in relation to the daily model output) which is halved recursively during the calculation in case of fraction variables becoming negative.

2.3.3. Atmospheric and Riverine N Input

The applied atmospheric and riverine TN inputs are crucial for the results of the TBNT analysis. Therefore, Figures 3A,B show time series of the annual atmospheric and riverine TN input into the North Sea domain, respectively. In both panels, the left y axes refer to the absolute TN input (in kt N a⁻¹). The right y axis in Figure 3A refers to the contributions of the atmospheric input into different subregions of the North Sea domain, relative to the total input into the domain: southern North Sea (▴; area of 1.7·10⁵ km²), northern North Sea (▾; 3.8·10⁵ km²), English Channel (◂; 0.9·10⁵ km²) and Kattegat/Belts (▸; 0.5·10⁵ km²). The southern North Sea represents the North Sea region south of 56 °N until transport section 5 in the Strait of Dover (see Figure 2). The northern North Sea represents the part of the TBNT domain north of 56 °N, including the Skagerrak north of section 4. The English Channel is the region between the EC boundary and section 5, while the Kattegat/Belts region is limited by the BS boundary and section 4. The total area of the North Sea domain is 6.9·10⁵ km².

Figure 3A shows that the atmospheric input (solid line) decreases slowly throughout 2000–2014, with a relative reduction of about 13% relating to the observed trend (Schöpp et al., 2003). The atmospheric deposition yields 490 kt N a⁻¹ in 2000 and stays around 485 kt N a⁻¹ during 2001–2004. From 2005 on, it shows a steady decrease reaching 430 kt N a⁻¹ in 2014.

The relative contributions of the atmospheric TN input into the different subregions of the TBNT domain (dotted lines) are stable throughout 2000–2014. This relates to the fact that only the total, domain-wide N deposition decreases, while the spatial patterns do not change. The largest part of about 44.7% is deposited in the northern North Sea. About 31.6% enter the southern North Sea, whose area is less than half of that of the northern North Sea region, emphasizing the remarkable spatial differences in atmospheric N deposition. The relative contributions to the English Channel region and the Kattegat/Belts are about 13.6% and 10.1%, respectively.

The absolute annual riverine TN input into the TBNT domain (Figure 3B; solid line) ranges between 775 kt N a⁻¹ in 2005 and 1266 kt N a⁻¹ in 2002. An apparent step change is seen between 2002 and 2003. The relative contributions of the different river groups to riverine TN input (dotted lines) show that the NL-1 rivers (incl. Rhine and Meuse; see Table 1) are the largest TN source during most years with values between 18.1% and 26.3%. The UK-2 rivers, located at the British east coast, are of similar importance ranging between 17.3% and 25.4%. The DE rivers constitute the third largest source during most years with contributions of 9.9–22.1%. The FR rivers contribute 10.6–18.2%. The UK-1 and other rivers show values between 6.1–13.1% and 5.7–14.1%, respectively, with higher values in the period after 2008 compared to the earlier period. This indicates that the relative importance of the different rivers on the North Sea system changed during this period. The remaining rivers (NL-2, BE, and NO) are less important (<5%).

The relative contributions show that the drop in TN input from 2002 to 2003 mainly relates to decreasing inputs from the DE and NL-1 rivers. Interestingly, the DE rivers show particularly high relative contributions in 2002 and 2010 (21.3% and 22.1%, respectively). Both relate to Elbe flood events in August 2002 (e.g., Ulbrich et al., 2003) and late summer 2010 (e.g., Kienzler et al., 2015; Philipp et al., 2015), respectively.

It should be noted that average annual PON:TN ratios for most groups are in the range of 0.16–0.19. Only the UK rivers show lower ratios around 0.08 and the NL-2 and NO rivers show comparably high ratios of about 0.45 and 0.47, respectively, though corresponding to small overall N loads. The standard deviation of the annual ratios is <0.05 for all groups. Thus, inorganic N inputs dominate the ecosystem dynamics.

3. Results

This section first provides a brief validation of simulated volume transports as well as spatial distributions of TN from different sources and minimum O₂ concentrations, simulated by the applied model and TBNT framework. By this, we demonstrate that the applied setup reproduces the main features of the North Sea ecosystem relevant to this study. This yields the basis for our analysis of the influence of N inputs on the O₂ dynamics.

3.1. Validation

3.1.1. Volume Transports

Table 2 presents volume transports across defined sections simulated by HAMSOM and corresponding transports from other modeling studies (Ozer, 2011; Winther and Johannessen, 2006) and observations (Winther and Johannessen, 2006). The 11 sections used for the comparison are indicated in Figure 2 (arrows indicate direction of positive net flow/inflow). They were selected with respect to the North Sea general circulation and data availability. Ozer (2011) conducted a model intercomparison on the North Sea circulation in 2008. Thus, the corresponding HAMSOM transports across these sections are calculated for 2008 only. For the other sections, given values represent statistics over the 15-year period 2000–2014. The provided observation data were compiled by Winther and Johannessen (2006) and originate from Otto et al. (1990), Rodhe (1996), Rydberg et al. (1996), and Danielssen et al. (1997).

HAMSOM reproduces the gross inflow from the North Atlantic across the northern shelf edge (Table 2, sections 1 to 3). Average gross inflow through the three sections simulated by HAMSOM adds up to 2.25 Sv. This is comparable to the results of Winther and Johannessen (2006) yielding 2.22 Sv and close to the upper end of the observed range (1.60–2.01 Sv). The ratio between the inflows through sections 2 and 3 is opposite to that for the transports simulated by Winther and Johannessen (2006) and to the literature values. However, it should be noted that these transports are highly sensitive to the position of the separation between sections 2 and 3. The higher minimum–maximum range in the HAMSOM transports (compared to Winther and Johannessen, 2006) likely relates to the longer analysis period. Regarding the outflow in the Norwegian Trench area (section 3), HAMSOM shows good agreement with the other studies.

The net transport through section 4 (Kattegat) corresponds with the results of Winther and Johannessen (2006) and observations. The same applies to the gross inflow through the English Channel (section 5), for which the average net inflow in 2008 matches with Ozer (2011). The comparison of the transports across the northern shelf edge (sections 1–3) with those through sections 4 and 5 emphasizes that the Atlantic inflow from the North constitutes the dominant inflow into the North Sea.

For sections 6 and 7, the net transports are in the range reported by Ozer (2011). Though, transports for section 6 are rather at its upper end. The negative net transport through section 6 implies southward flow in the southwestern North Sea driven by the Atlantic inflow from the North and freshwater discharge from the rivers at the British east coast. In contrast, the northwestward flow in the eastern part (section 7) is governed by the continental coastal current driven by the English Channel inflow and riverine freshwater especially from the large Dutch rivers.

For the net flow into and out of the German Bight (sections 8 and 9, respectively), HAMSOM yields values of about 0.09 Sv, slightly below the range provided by Ozer (2011). This likely relates to the bathymetry in combination with the vertical resolution resulting in a distinct step in the model bathymetry causing a northward turn of parts of the coastal current before crossing section 8 (not shown).

The negative net flow across section 10 in the western central North Sea simulated by HAMSOM indicates the mainly southward flow in this region. The average value of −0.44 Sv is in the range found by Ozer (2011). In the eastern central North Sea (section 11), the positive net flow derived from HAMSOM represents the generally northward flow in this region relating to the cyclonic recirculation of most of the southward flowing water masses north of Dogger Bank (Lenhart and Pohlmann, 1997). The transport is slightly above the range reported by Ozer (2011), which may relate to a lower inflow from the Baltic Sea indicated by the differences in the transports across section 4.

In summary, the HAMSOM volume transports across the outer boundaries and in the interior of the North Sea show good agreement with other models and observations and represent well the general cyclonic circulation (e.g., Otto et al., 1990). Therefore, the current fields simulated by HAMSOM provide a reliable basis for the TBNT analysis.

3.1.2. Spatial Distributions of TN from Different Sources

In order to give insight in the dispersal of N from different input sources in the North Sea, Figure 4 shows the mass-weighted average relative contributions of selected input sources to TN (i.e., the sum of all N-related model state variables incl. sediment-N) during 2000–2014: (Figure 4A) atmosphere, (Figure 4B) German rivers (DE), (Figure 4C) 1st group of Dutch rivers (NL-1), and (Figure 4D) 2nd group of British rivers (UK-2). By aggregating pelagic and benthic state variables the internal processes (e.g., uptake of $\text{N}{\text{O}}_3^{-}$ by diatoms) cancel out and the spatial distributions of TN from the different sources basically reflect the effect of lateral exchange processes. Only the loss of molecular N₂ due to benthic denitrification, and atmospheric N deposition have an additional influence. The individual N sources were selected in relation to existing TBNT studies on TN (OSPAR, 2010; Painting et al., 2013; Troost et al., 2013).

The relative contributions show that the southeastern North Sea is strongly affected by atmospheric N deposition (Figure 4A) as well as the selected riverine N sources (Figures 4B–D). In contrast, these input sources are of minor importance (<2.5%) in the entire northern and large parts of the central North Sea. In these regions, TN is dominated by the NA influence (not shown) which decreases from 100% in the off-shore regions of the northern North Sea to about 50% along the 50-m-isobath.

In the coastal northwestern North Sea, the UK-2 (see Figure 4D) and other rivers (not shown) have an increased influence. In the northeastern North Sea, the atmosphere (see Figure 4C) and the Norwegian rivers (NO; not shown) constitute relevant contributions besides the NA contribution.

The atmospheric contribution (Figure 4A) plays a relevant role in the southern central and southern North Sea, ranging between 5% in the southwestern North Sea and 19% near the Danish west coast. Furthermore, the contributions of the DE, NL-1, and UK-2 rivers (Figures 4B–D, respectively) to TN are high in the southern central and southern North Sea. In the vicinity of the river mouths of Rhine/Meuse (NL-1), Humber and Wash (UK-2) as well as Elbe and Weser (DE), the contributions of the corresponding rivers can reach values of almost 100%. However, their contributions decrease to <40–60% within a distance of about 100 km from each river mouth.

Though, the contribution of the NL-1 and UK-2 rivers (Figures 4C,D, respectively) extends far into the off-shore regions of the southern North Sea. The former reaches values well above 10% in almost the entire southeastern North Sea, while the latter stays well above 10% in the entire southwestern part. The DE inputs (Figure 4B) extend northward along the German and Danish west coast. This illustrates the cyclonic general circulation of the North Sea.

The relevant contributions of the selected TN sources are accompanied by a steady decrease in the NA and EC contributions (not shown) from about 60% in the off-shore parts to 5–10% in the coastal regions between the mouths of Rhine/Meuse and the Danish northwest coast.

3.1.3. Minimum Bottom O₂ Concentrations

Große et al. (2016) showed that the applied model setup reproduces the main spatial and temporal features of North Sea O₂, only slightly overestimating near-bottom O₂ concentrations. Therefore, Figure 5 only shows the spatial distribution of the overall minimum simulated bottom O₂ concentration during 2000–2014. By this, a first insight in the concurrence of N from individual input sources and low O₂ concentrations is provided.

Lowest simulated O₂ concentrations occur in the region between 54.5–56.5 °N and 4–7.5 °E. This is in good agreement with recurrent observations of O₂ deficiency in this region in recent years (Figure 14 in Topcu and Brockmann, 2015). Minimum O₂ concentrations of 5.2 mg O₂ L⁻¹ are found in region 1 (see Figure 5) which, therefore, is used for a detailed analysis of the O₂ dynamics in relation to N inputs. Region 1 is hereafter referred to as “O₂ deficiency zone” (ODZ). The comparison with the spatial patterns of TN from individual sources (see Figure 4) indicates that this region is under strong influence of the DE rivers, but also the NL-1 and UK-2 rivers. The atmosphere and the NA also play a relevant role.

Low simulated O₂ concentrations also occur directly north of Dogger Bank which matches with observations (Weston et al., 2008; Greenwood et al., 2010; Queste et al., 2016). The same applies to the Oyster Grounds site (see Figure 5, region 2) where simulated minimum near-bottom concentrations yield 6–6.5 mg O₂ L⁻¹. As the Oyster Grounds are well known being susceptible to low O₂ conditions (de Wilde et al., 1984; Peeters et al., 1995; Weston et al., 2008; Greenwood et al., 2010), this region is also used for a detailed analysis. At this site, the comparison with Figure 4 indicates a strong influence of the NL-1 and UK-2 rivers, and the NA (not shown). The atmospheric contribution is lower than at the ODZ site. It should be noted that the term “Oyster Grounds” used in this study refers explicitly to the defined region 2 which was chosen in relation to the study site of Greenwood et al. (2010).

3.2. O₂ Consumption Related to N from Individual Sources

Based on the good agreement of the presented results with existing observations and modeling studies, the following analyses focus on the influence of N from selected input sources on near-bottom GOC.

3.2.1. Spatial Distributions

Figure 6A shows the spatial distributions of annual GOC in the model bottom layer averaged over 2000–2014. The corresponding average relative contributions related to selected N input sources are given in Figures 6B–F: (Figure 6B) boundaries (NA, EC, and BS), (Figure 6C) atmosphere, (Figure 6D) DE rivers, (Figure 6E) NL-1 rivers, and (Figure 6F) UK-2 rivers.

GOC shows a general South-North gradient with high values (>0.3 kg O₂ m⁻² a⁻¹) in the coastal southern North Sea and low values (<0.1 kg O₂ m⁻² a⁻¹) in the deep northern parts. The lower values in the deep North Sea regions relate to lower near-surface NPP and greater bottom depths, resulting in less organic matter reaching the bottom. Highest GOC occurs in the vicinity of large river mouths (e.g., Elbe) due to riverine input of organic matter additionally enhancing local GOC. Here, values of up to about 1.3 kg O₂ m⁻² a⁻¹ are simulated.

It can be seen that the coastal regions of highest GOC do not correspond to the regions of lowest bottom O₂ concentrations (see Figure 5). This relates to the strong tidal mixing inhibiting long-lasting seasonal stratification which constitutes the prerequisite for the evolution of low O₂ conditions (Greenwood et al., 2010; Große et al., 2016). In the off-shore regions of the southeastern North Sea average annual GOC ranges between 0.1 kg O₂ m⁻² a⁻¹ and 0.2 kg O₂ m⁻² a⁻¹. Here, GOC increases with decreasing water depth due to higher organic matter availability near the seafloor (Große et al., 2016).

The spatial distributions of source-specific relative contributions to GOC (see Figures 6B–F) basically represent the patterns shown for TN (see Figure 4). Thus, the input of N from the boundaries (see Figure 6B), especially the NA, dominates most off-shore regions of the North Sea, while the atmosphere and the riverine contributions from NL-1 and UK-2 (see Figures 6C,E,F, respectively) play an important role for GOC in the southern North Sea. The DE rivers' influence on GOC is confined to the German Bight and Danish west coast (see Figure 6D).

Minor differences can be seen with respect to the relative importance of the atmosphere which shows a slightly higher contribution to GOC in most off-shore areas compared to its contribution to TN. This likely relates to the stronger influence of atmospheric N deposition on NPP as deposited inorganic N is immediately available to primary producers. At river mouth locations, the corresponding river groups show slightly higher contributions to GOC than to TN which underlines the relevance of riverine organic matter for local GOC.

3.2.2. Inter-Annual Variability in the ODZ

Figure 7 provides a more detailed insight in the influence of N from individual sources on GOC in the region most susceptible to O₂ deficiency. It shows the time series of average daily, near-bottom GOC at the ODZ site (see Figure 5, region 1) during the stratified period (t_strat) of the individual years 2000–2014, and average values over t_strat of the entire period. t_strat is derived from simulated T following Große et al. (2016). The entire period from the very first day of stratification until the very last day (incl. intermittences) is used for averaging. Consequently, average daily GOC provides a comparable measure for analyzing inter-annual variations in the influence of N from individual sources on the O₂ dynamics, independent of the duration of seasonal stratification. For a more quantitative description, Table 3 (left-hand side) provides the corresponding average (and standard deviation) values for t_strat, average daily GOC and the relative contributions of individual N sources to average daily GOC for the entire period 2000–2014, and for the years 2002 and 2010. All values refer to a bottom layer of 6.5 m average thickness.

Figure 7 and Table 3 show that 2002 represents an exceptional year with the highest daily GOC (0.61 g O₂ m⁻² day⁻¹). This high GOC coincides with the longest t_strat of 186 days, resulting in the lowest bottom O₂ concentrations in the vicinity of the ODZ site in 2002 (Große et al., 2016). In contrast, 2010 constitutes the year of lowest GOC (0.31 g O₂ m⁻² day⁻¹) resulting in higher bottom O₂ concentrations. Interestingly, 2011 also shows high daily GOC of 0.56 g O₂ m⁻² day⁻¹, though, O₂ deficiency does not occur in 2011 as a result of comparably short t_strat of only 135 days (not shown). This puts emphasis on the crucial role of stratification for the evolution of O₂ deficiency in the off-shore North Sea.

The NA contribution constitutes the major influence on N-related GOC at the ODZ site, with an average contribution of 39.8% to daily GOC and maximum contributions of up to 54.2%. This maximum is reached in 2000, when daily GOC reaches its second highest value of 0.58 g O₂ m⁻² day⁻¹. However, it can be seen that years of high contributions by the NL-1 or DE rivers can also show very high GOC (2002 and 2011, respectively). The NL-1 contribution is highest in 2002 resulting in 15.4%, while the DE contribution is highest in 2011 (20.1%). This can be attributed to the Elbe flood in January 2011 (e.g., Kienzler et al., 2015; Mudersbach et al., 2017). Interestingly, the Elbe flood in August 2002 (Ulbrich et al., 2003; Kienzler et al., 2015; Mudersbach et al., 2017) is not reflected in a higher DE contribution in 2002, as the related loads reached the ODZ site after the growing season. However, it explains the large DE contribution of 18.2% in 2003 which is the second largest during the entire period.

The atmosphere constitutes the second largest contribution to GOC at the ODZ site, while the DE and NL-1 rivers constitute the major riverine sources yielding 11.0% and 10.0%, respectively (see Table 3). The British rivers (UK-1 and UK-2) also play a relevant role, adding up to 9.2%. The other sources are of minor importance with average contributions of <5%. However, they can show higher values during individual years (see maximum values in Table 3). The minimum–maximum ranges further indicate that especially the DE and NL-1 contributions as well as the NA contribution can show high variations between individual years. The average integrated riverine contribution to GOC of 39.1% is in the same order as the NA contribution and only about 5% less than the integrated contribution from the boundaries (NA, EC, and BS). This illustrates the great importance of riverine N inputs on O₂ deficiency at the ODZ site, especially, from the German, Dutch, and British rivers.

3.2.3. Inter-Annual Variability in the Oyster Grounds

Figure 8 shows the time series of average daily, near-bottom GOC at the Oyster Grounds site (see Figure 5, region 2) analogous to Figure 7. The corresponding values are given in Table 3 (right-hand side).

Average daily GOC is generally lower than at the ODZ site, which only partly relates to the smaller bottom layer thickness of 5.5 m. As for the ODZ site, minimum and maximum GOC rates occur in the years 2010 (0.23 g O₂ m⁻² day⁻¹) and 2002 (0.39 g O₂ m⁻² day⁻¹), respectively. However, some differences between the two sites can be seen for the other years, e.g., in 2011, which constitutes the year of only seventh strongest GOC in the Oyster Grounds.

The contributions of the individual N sources to GOC reveal the causes for these regional differences. First, the DE contribution is very small in most years as the DE rivers are downstream the Oyster Grounds relative to the general circulation. Consequently, the high N inputs from these rivers in 2011 do not affect GOC in the Oyster Grounds region. Second, the absolute NA contribution, which is the largest among all contributions, is clearly lower (average of 0.13 g O₂ m⁻² day⁻¹) than at the ODZ site (average of 0.19 g O₂ m⁻² day⁻¹). However, the relative NA contribution is about the same as at the ODZ site (39.3%) due to the lower overall GOC. In addition, the absolute atmospheric contribution is lower in the Oyster Grounds (average of 0.04 g O₂ m⁻² day⁻¹ compared to 0.06 g O₂ m⁻² day⁻¹ at the ODZ site), which also results in a slightly lower relative contribution of 13.1%.

Table 3 reveals further differences for the remaining N sources compared to the ODZ site. The influence of the NL-1 and UK-2 rivers is visibly higher, yielding average relative contributions of 14.6% and 12.9%, respectively. These values correspond to a similar average absolute NL-1 contribution (0.05 g O₂ m⁻² day⁻¹) as for the ODZ site and an increased average absolute contribution of the UK-2 rivers (0.04 g O₂ m⁻² day⁻¹ compared to 0.03 g O₂ m⁻² day⁻¹). Consequently, years of high contributions of either of the two N sources result in high GOC in the Oyster Grounds (see Figure 8; e.g., 2002, 2003, 2013). The remaining contributions are less important for GOC in the Oyster Grounds. Though, minor increases in the relative contributions can be seen in basically all contributions, except for the NL-2 and NO rivers, and the BS. This relates to the downstream locations of these sources relative to the Oyster Grounds.

The average integrated riverine contribution to N-related GOC in the Oyster Grounds adds up to 40.9%, with 30% attributed to the Dutch and British rivers. This is even higher than the NA contribution and only about 5% less than the integrated contribution from the boundaries, putting emphasis on the importance of riverine N on the O₂ dynamics in the Oyster Grounds.

4. Discussion

This study, for the first time, presents a detailed quantitative analysis of the O₂ dynamics of the North Sea in relation to N inputs from various external sources, including rivers, the atmosphere and adjacent seas. By this, the presented approach constitutes a valuable tool for the source-oriented approach advocated by OSPAR (1999).

There have been earlier studies on the relationship between river nutrient loads and O₂ deficiency (e.g., Brockmann et al., 1988; Topcu and Brockmann, 2015). However, these studies rely on statistical analyses and the qualitative description of the interactions controlling the evolution of O₂ deficiency. Other studies investigated the effect of river nutrients on the North Sea biochemistry based on nutrient reduction scenarios (e.g., Skogen et al., 2004; Lenhart et al., 2010). Though, only Lenhart et al. (2010) considered implications for the O₂ conditions in the North Sea based on a model comparison.

However, these scenario-based approaches do not allow for the analysis of the current state of the ecosystem, that is represented by a reference scenario with realistic forcing data, with respect to the influence of nutrients from individual sources on the O₂ dynamics. Furthermore, the assessment of the influence of individual sources would require an individual scenario for each source, which is time-consuming and costly. In contrast, the present method requires only one model simulation and one TBNT analysis dealing with all selected N input sources at once.

The source-specific information on GOC obtained from the expanded TBNT method can be used for speeding up the iterative process of finding optimal N reduction targets using the linear optimization method described by Los et al. (2014). Nevertheless, actual N reduction scenarios will still be necessary in order to assess the actual impact of the derived reductions targets, due to the non-linear ecosystem interactions.

With respect to the presented results, the simulated volume transports show good agreement with other models and literature values. The same applies to the spatial distributions of TN from individual sources which correspond to results of other North Sea TBNT studies (OSPAR, 2010; Painting et al., 2013; Troost et al., 2013), giving confidence for the applicability of a TBNT post-processing software. Consequently, the analysis of GOC related to N from these sources can be assumed to yield valid information. However, the applied model setup has some limitations which require discussion.

4.1. Limitations of the Model Forcing and Setup

The forcing provided for the different N input sources (rivers and atmosphere) has a strong influence on the results of the TBNT analysis. Although the river data is compiled to the best knowledge and represents the state of the art of North Sea river data, it is based on assumptions as a result of limited data availability. An upscaling based on the distance between gauge station and river mouth and catchment area size is applied to the British rivers (van Leeuwen, pers. comm.). Here, the uncertainty can become high in case of a large distance. In addition, for some rivers only monthly data exist which have to be interpolated to daily time series, partly with support of climatologies which may not fully represent the recent state. Furthermore, estuarine retention of N and P is not considered, although they were found to be non-negligible, e.g., in the Elbe estuary (Behrendt, 1996; Dähnke et al., 2008).

Furthermore, point sources (e.g., sewage treatment plants) in the estuaries — downstream the nutrient monitoring locations used for the load calculation (Pätsch et al., 2016) — are not considered in the river data, except for the British rivers. For the major Dutch rivers (excl. the smaller North Sea Canal) nutrient concentrations are measured right at the floodgates, i.e., point sources are implicitly considered. For the German rivers, which together with the former two are the most important riverine N sources, only limited information on point sources are available. However, Dähnke et al. (2008) found no evidence of a significant influence of point sources on $\text{N}{\text{O}}_3^{-}$ concentrations in the Elbe estuary, based on stable isotope analyses. A catchment modeling study for the Weser shows that only about 0.1 kt N a⁻¹ enter the surface waters along the Weser estuary (Heidecke et al., 2015). Hence, it can be assumed that neglecting point sources in the estuaries is a reasonable constraint.

For the atmospheric N deposition the data availability also constitutes a relevant constraint. In order to run a realistic simulation, the data needs to be provided for the entire simulation period and model domain. This can be obtained by temporal extrapolation and spatial interpolation to the model grid, based on the original N deposition data from EMEP. However, this original data as well as the spatio-temporal processing are subject to uncertainties. In this study, annual EMEP data was used, though, also monthly data are available (Desmit et al., 2015). As atmospheric N deposition in the North Sea region experiences a distinct annual cycle with higher summer and lower winter values (e.g., Neumann et al., 2016), the use of these monthly data could be worthwhile for future analyses. Differences in the applied atmospheric deposition data may also partly explain the minor differences between the results for the atmospheric contribution to TN of this study and those of Troost et al. (2013), who found a higher atmospheric contribution in the southern North Sea.

Another N-related aspect not considered in ECOHAM is the exchange between the North Sea and the Wadden Sea. Postma (1981) describe the Wadden Sea as a relevant sink for organic and source of inorganic matter. Other studies reported the Wadden Sea constitutes an important sink of N due to benthic denitrification (e.g., Kieskamp et al., 1991; van Beusekom and de Jonge, 1998; van Beusekom et al., 1999). Both would affect N availability in the open North Sea, especially, upstream of the ODZ site. Therefore, the Wadden Sea should be included in future studies.

With respect to the TBNT results, we consider the NA contribution to TN (and GOC) in the southeastern North Sea to be slightly overestimated. We conclude this from a comparison of the source-specific TN contributions of the other most relevant sources (atmosphere, DE, NL-1, and UK-2) to a box in the inner German Bight (not shown) calculated in our study with other TBNT studies on TN (e.g., OSPAR, 2010; Painting et al., 2013; Troost et al., 2013). This may relate to the limited spatial resolution of about 20 km that does not allow for the representation of sub-mesoscale features such as eddies (e.g., Johannessen et al., 1989; Røed and Fossum, 2004; Albretsen, 2007).

4.2. Limitations of the Model Intrinsics and the TBNT Method

With respect to the hydrodynamics, the component-upstream advection scheme applied by HAMSOM could be an additional cause for the overestimation of the NA contribution as it results in high numerical diffusion (Lenhart and Pohlmann, 1997). However, the identification of the actual cause goes beyond the scope of this study.

Regarding the model biogeochemistry, it should be noted that the benthic re-mineralization tends to be slightly underestimated by ECOHAM (Große et al., 2016). On the one hand, this can affect the overall GOC. On the other hand, the contributions of individual N sources to GOC may be affected as the influence of sources with high contributions to benthic re-mineralization would be amplified relative to sources with only low contributions. In addition, the rather simple sediment module does not allow for long-term storage or burial of organic matter in the sediment. However, this matches with observational studies that found burial of N and C in the North Sea sediment being only small relative to annual cycles of NPP and respiration (Wirth and Wiesner, 1988; de Haas et al., 1997; Brion et al., 2004).

In relation to the benthic processes, it should also be noted that simulated daily benthic denitrification rates in the inner German Bight are up to 4 times higher than those recently reported by Marchant et al. (2016). However, this only applies to few very near-shore areas in the model. In the off-shore regions nitrification rates range between 0.4 and 2 mmol N m⁻² day⁻¹, which is in the range they found. Hence, we rate the additional N loss in our simulation as acceptable as it only occurs in a very small region. In addition, the here applied simple benthic module only includes the effect of O₂ availability in the pelagic bottom layer on benthic processes. As anoxic conditions do not occur in the pelagic North Sea, this has no effect on the benthic processes. Thus, the incorporation of a more complex benthic module, that resolves vertical O₂ gradients within the sediment, would be worthwhile.

Furthermore, ECOHAM applies constant C:N:P ratios to phytoplankton, zooplankton and bacteria (Lorkowski et al., 2012; Große et al., 2016). Though, it is known that elemental ratios vary between different phytoplankton species (e.g., Geider and La Roche, 2002; Quigg et al., 2003; Malzahn et al., 2010) as well as under varying environmental conditions (e.g., Rhee, 1978; Goldman et al., 1979). In fact, different elemental ratios, e.g., in phytoplankton, would impact on the entire ecosystem including organic matter production and, thus, on the O₂ dynamics. For instance, unpublished sensitivity studies showed that organic matter export below 30 m varies by up to 13% with varying stoichiometry (Redfield ± 50%). Consequently, this may slightly change the relative contributions of the different N sources to GOC in relation to temporal variability in these contributions combined with temporal variability in the elemental ratios. However, in order to assess the impact of different stoichiometry a comprehensive sensitivity analysis would be required which is beyond the scope of this study.

With respect to the expansion of TBNT by the link to the O₂ cycle, the use of the mediators in case of C-based O₂-affecting processes (e.g., NPP) constitutes a key element. Therefore, discussion of the constraints of this approach is required. The actual values of the C-based O₂-affecting processes are proportional to the concentration of the C-based mediator state variable (see Appendix A.1–A.4). The use of the mediator implies that the concentration of the C-based mediator variable is proportional to the corresponding N-based mediator variable. This evidently can be applied in the case of a fixed C:N ratio as used for phytoplankton, zooplankton and bacteria in the present study.

Conceptually, this also holds for a varying stoichiometry as used for benthic detrital matter (i.e., sediment-N/-C) in ECOHAM. Higher concentrations of sediment-N allow for higher concentrations of sediment-C and, thus, benthic O₂ consumption. However, temporal variations in the C:N ratio in benthic detrital matter in combination with temporal variations in the contributions of individual N sources to sediment-N can affect the contributions of these sources to benthic O₂ consumption and, thus, to near-bottom GOC.

A quantitative analysis of this effect would require a comprehensive sensitivity study on the effect of different C:N ratios on the TBNT results. Here, we will only discuss effects of the mediator use in relation to potential event-driven changes in the C:N ratio in benthic detrital matter.

In the off-shore regions of the southern and southeastern North Sea — which are in the focus of this study — benthic organic matter originates mainly from local NPP (Große et al., 2016). Thus, organic matter input from rivers especially during the summer period does not need to be considered. Davies and Payne (1984) observed C:N ratios of 9–18 in particulate organic matter 2 m above the sea floor in the northern North Sea during the spring bloom 1976. For the southern North Sea off Belgium, Tungaraza et al. (2003) report C:N ratios of 5.7–11.2 in particulate organic matter, depending on the phytoplankton species composition. This indicates that even over short periods of time (days to weeks) C:N ratios in organic matter may vary by a factor of 2. If such short-term variations coincide with strong variations in the contributions by individual N sources to detrital matter, this may result in an under-/overestimation of the influence of individual N sources on GOC. Analogously, such effects may occur, e.g., during events of strong rainfall over sea during the growing season, which enhance wet N deposition and thereby enable new production due to increased $\text{N}{\text{O}}_3^{-}$ availability (e.g., Dugdale and Goering, 1967).

Regarding the here presented analysis of source-specific GOC at the ODZ and the Oyster Grounds sites it should be noted that the C:N ratios in benthic detrital matter show only very minor variations, ranging from 7.5 to 7.9, during the growing season. As this period basically corresponds to the summer stratified period, which is in the focus of this study, the effect of changing C:N ratios on the source-specific contributions to GOC is considered to be small. However, considering model applications using variable phytoplankton stoichiometry additional efforts would be required for further improving the mediator approach. Here, it may be feasible to use the actual C:N ratios and relative contributions by the individual N sources in both, detrital matter that was previously sedimented and newly sedimented matter, as a weighing factor. By this potential errors in the current approach could be minimized.

Besides the effects of varying C:N ratios, differences in the re-mineralization rates of C (0.028 day⁻¹) and N (0.0333 day⁻¹) may affect the results of the expanded TBNT due to a temporal delay between C and N re-mineralization (and benthic denitrification). However, in the present study the time series of sediment-N and -C (not shown) reveal maximum temporal delays of only 2–3 days. Therefore, we consider this to have only a minor impact on the results.

In relation to the mediators, it has to be stated clearly that the results of this study can only be interpreted with respect to N, i.e., they do not reflect the overall influence of the selected input sources on the O₂ dynamics. Labeling C or P from individual input sources will likely result in different relative contributions. On the one hand, this may relate to changes in the inputs by the different sources relative to each other (e.g., different ratios between P loads of individual rivers compared to N loads). On the other hand, selected processes (e.g., nitrification) cannot be considered in case of P labeling, while other processes such as P release from the sediment under anoxic conditions would play an additional role when considering P (e.g., Hupfer and Lewandowski, 2008). However, this process is not included in ECOHAM as it is designed primarily to represent the O₂ cycle in relation to N and C. As N constitutes the limiting nutrient among these two with respect to organic matter production, we focused our analysis on N. In order to assess the overall impact of nutrients (i.e., N and P) from individual sources on the North Sea O₂ dynamics, an analogous analysis for P would be required. However, this is beyond the scope of this study. In relation to that, it also has to be noted that the TBNT analysis does not replace actual nutrient reduction scenarios, but rather constitutes a useful tool for speeding up the iterative process of finding an optimal nutrient reduction setup (Los et al., 2014).

Furthermore, the presented method cannot be applied to formulations of benthic O₂ consumption only based on the O₂ concentration itself, e.g., that by Hetland and DiMarco (2008), as those cannot be linked to external nutrients.

Due to these uncertainties and constraints, we highly recommend applying the presented method to different models and model setups. This would provide a range of the contributions of the individual N sources on the O₂ dynamics of the North Sea, giving a more reliable estimate of their actual influence. The only requirement for such approach is the process-based parametrization of the O₂-affecting processes in that model, such that a mediator substance provides a link between the labeled element and the O₂ cycle.

5. Conclusions

The here presented expansion of the TBNT method and the corresponding analysis of source-specific GOC represent the first approach for quantifying the influence of N inputs from different sources on the O₂ dynamics of a marine ecosystem. Our study focused on the North Sea and our analysis constitutes a major advancement toward the source-oriented approach for the ecological management of O₂ deficiency in the North Sea (OSPAR, 1999). Besides this, the described method can be applied analogously to different elements (e.g., P or C), ecosystems and models. Therefore, it could be used for the investigation of other ecosystems affected by eutrophication-related O₂ deficiency, e.g., the Baltic Sea (e.g., Rosenberg, 1985) or the northern Gulf of Mexico (e.g., Rabalais and Turner, 2001).

The present analysis focused on the general spatial patterns and inter-annual variability in the regions most susceptible to O₂ deficiency. Though, the applied method can also be used for analyzing the influences of individual sources of N (and other nutrients) on different time scales (days to decades).

While GOC in the northern and northern central North Sea is dominated by the N inputs from the North Atlantic, GOC in the southern central and southern North Sea is under the combined influence of the North Atlantic (NA), the English Channel (EC), the atmosphere and the large riverine N sources. Among the riverine N sources, especially, the large Dutch rivers (NL-1; incl. Rhine and Meuse) and the rivers at the British east coast (UK-2) have a relevant influence on GOC in the off-shore regions of the southern North Sea, while the German rivers dominate the southeastern coastal North Sea.

The analysis of GOC at the two study sites confirmed that seasonal stratification constitutes the prerequisite for the evolution of O₂ deficiency in the off-shore regions of the southern North Sea (Große et al., 2016). It further revealed that riverine N inputs play a major role for the bottom O₂ dynamics, with average relative contributions of about 40% during 2000–2014. Depending on their timing, flood events in the German or Dutch rivers can increase the riverine contribution to even more than 50%. In addition, atmospheric N contributes about 15% to the overall GOC at both sites. As riverine and atmospheric inputs are strongly anthropogenically affected (Paerl, 1997; Schöpp et al., 2003) this emphasizes the large human impact on the O₂ conditions in the North Sea regions most susceptible to O₂ deficiency.

Consequently, thoroughly defined riverine N reductions, e.g., under the European Union's Water Framework Directive (EU, 2000) likely constitute a potent means for improving the O₂ conditions in the southern North Sea. In combination with potential future reductions in atmospheric N deposition — as adopted within the Gothenburg Protocol (UNESC, 1999) — such reductions may improve the North Sea O₂ conditions or mitigate potential negative effects of climate change.

In this context it is important to note that also future changes in the physical environment, e.g., in stratification or seawater T, will have a relevant impact on the O₂ conditions of the North Sea as they play a key role in the O₂ dynamics (Greenwood et al., 2010; Große et al., 2016). However, our results emphasize that future scenarios for the North Sea require realistic estimates for both, the climatic conditions and the riverine and atmospheric N loads in order to provide a likely picture of the future North Sea O₂ conditions.

Author Contributions

The original manuscript was conceptualized and set up by FG and H-JL. The simulations were set up and conducted jointly by FG, MK, and JP. All analyses were conducted by FG. All authors contributed equally to the reworking of the original manuscript toward the submitted version.

Funding

This study received funding from the German Environmental Protection Agency (UBA), in the frame of the project “Implementation of Descriptor 5 Eutrophication to the “MSFD”, SN: 3713225221, and the Cluster of Excellence “CliSAP” (EXC177), Universität Hamburg, funded by the German Science Foundation (DFG).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the three reviewers for their helpful comments and constructive criticism on the original manuscript. We would like to thank Sonja van Leeuwen from Cefas for providing freshwater discharge and nutrient load data for the major rivers across Europe. We further thank Jerzy Bartnicki from EMEP for providing atmospheric N deposition data. The model simulations were conducted on “Mistral,” the Atos bullx DLC B700 mainframe at the German Climate Computing Center (DKRZ) in Hamburg. The methods used in this study were published as part of the dissertation by Große (2017).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2017.00383/full#supplementary-material

References