Implications of different nitrogen input sources for potential production and carbon flux estimates in the coastal Gulf of Mexico (GOM) and Korean Peninsula coastal waters

The coastal Gulf of Mexico (GOM) and coastal sea off the Korean Peninsula (CSK) both suffer from human-induced eutrophication. We used a nitrogen (N) mass balance model in two different regions with different nitrogen input sources to estimate organic carbon fluxes and predict future carbon fluxes under different model scenarios. The coastal GOM receives nitrogen predominantly from the Mississippi and Atchafalaya rivers and atmospheric nitrogen deposition is only a minor component in this region. In the CSK, groundwater and atmospheric nitrogen deposition are more important controlling factors. Our model includes the fluxes of nitrogen to the ocean from the atmosphere, groundwater and rivers, based on observational and literature data, and identifies three zones (brown, green and blue waters) in the coastal GOM and CSK with different productivity and carbon fluxes. Based on our model results, the potential primary production rate in the inner (brown water) zone are over 2 gC m⁻² d⁻¹ (GOM) and 1.5 gC m⁻² d⁻¹ (CSK). In the middle (green water) zone, potential production is from 0.1 to 2 (GOM) and 0.3 to 1.5 gC m⁻² d⁻¹ (CSK). In the offshore (blue water) zone, productivity is less than 0.1 (GOM) and 0.3 (CSK) gC m⁻² d⁻¹. Through our model scenario results, overall oxygen demand in the GOM will increase approximately 21 % if we fail to reduce riverine N input, likely increasing considerably the area affected by hypoxia. Comparing the results from the USA with those from the Korean Peninsula shows the importance of considering both riverine and atmospheric inputs of nitrogen. This has direct implications for investigating how changes in energy technologies can lead to changes in the production of various atmospheric contaminants that affect air quality, climate and the health of local populations.

1 Introduction

Industrial expansion and anthropogenic emissions are major factors leading to increased coastal productivity and potential eutrophication (Sigman and Hain, 2012). Coastal primary production is controlled largely by nitrogen (N) and phosphorus (P), and the relative supply of each determines which element limits production (Paerl, 2009); freshwater inputs and the distance from sources such as river mouths are also important (Dodds and Smith, 2016). Changes in nutrient loading from airborne, river-borne and groundwater sources can also affect which element limits coastal productivity (Sigman and Hain, 2012). Most coastal regions are N-limited; however, at certain times conditions can change from N-limited to P-limited (Dodds and Smith, 2016; Howarth and Marino, 2006). Sylvan et al. (2006), for example, suggested that the coastal Gulf of Mexico (GOM), especially near the Mississippi River delta mouth, is P-limited at certain times.

Several studies have shown that increasing atmospheric nitrogen deposition (AN-D) is contributing to ocean production globally, including to eutrophication, and is potentially of future importance in the GOM (Cornell et al., 1995; Doney et al., 2007; Duce et al., 2008; He et al., 2010; Kanakidou et al., 2016; Kim, 2018; T. W. Kim et al., 2011; Lawrence et al., 2000; Paerl et al., 2002). Recently, T. W. Kim et al. (2011), using a model simulation, showed that AN-D controls approximately 52 % of the coastal productivity in the Yellow Sea. Global NO_x emissions have increased but appear to be changing differently in the USA and Asia (J. Y. Kim et al., 2010; Luo et al., 2014; Shou et al., 2018; Zhao et al., 2015), and may affect not only coastal productivity but also global total nitrogen budgets. This study uses a box model to define potential carbon fluxes based on different nitrogen input sources in two different regions, the coastal Gulf of Mexico (GOM) and the coastal sea off the Korean Peninsula (CSK). The GOM and CSK were selected in this study because while the major input source to the coastal ocean in both regions is riverine, the AN-D and submarine groundwater discharge (SGD) are considerably more important in the CSK region (Wade and Sweet, 2008; Zhao et al., 2015).

Most previous model studies in the GOM have been used to predict the size of the hypoxic zone (e.g., Fennel et al., 2006, 2011, 2013; Green et al., 2008; Hetland and DiMarco, 2008; Justic et al., 2002; Scavia et al., 2004; Turner et al., 2006, 2008), although Bierman et al. (1994) used a mass balance model to estimate carbon flux and oxygen exchange. The mass balance model is a useful tool to calculate nutrient or carbon fluxes and to estimate production in the coastal ocean (J. S. Kim et al, 2010; G. Kim et al., 2011), and such models have been successfully used in many regions and individual coastal systems to estimate ecosystem metabolism, e.g., in the Patuxent River estuary of the Chesapeake Bay (Hagy et al., 2000; Testa et al., 2008) and in the LOICZ (Land Ocean Interactions in the Coastal Zone) project (e.g., Ramesh et al., 2015). However, there are few such model studies in the GOM and CSK. All previous models for the GOM and the CSK have considered only riverine N as the predominant input source, and AN-D as an input in either region has not been considered.

In this study, we aimed to (1) build a mass balance model considering not only riverine N input but also airborne and groundwater-borne N; (2) use it to calculate potential primary production in the three regions defined by Rowe and Chapman (2002, henceforth RC02, see next section) and their associated coastal productivity; and (3) use the mass balance model to test the RC02 hypothesis. Because RC02 did not quantify their model with nutrient data, and because this model has not been applied to another region, we tested the RC02 hypothesis using data from both the GOM and the CSK that include low salinity samples. We used historical data from the mid-western part of the CSK and evaluated the theoretical model of RC02 in both areas where freshwater with high terrestrial nutrient input mixes into the coastal ocean.

2 Study areas

The Louisiana–Texas (LATEX) shelf in the northern Gulf of Mexico is affected by coastal nutrient loading, leading to hypoxia, coming from two major terrestrial sources (the Mississippi and Atchafalaya rivers that together form the Mississippi–Atchafalaya river system, MARS). These two major rivers have different nutrient concentrations. The Gulf of Mexico (GOM) is a semi-enclosed oligotrophic sea and the MARS is the major source of nutrients and freshwater to the northern GOM (Alexander et al., 2008; Rabalais et al., 2002; Robertson and Saad, 2014). The MARS drains 41 % of the contiguous United States (Milliman and Meade, 1983) and discharges approximately 20 000 m³ s⁻¹, or about 60 % of the total freshwater flow, (about 10.6×10¹¹ m³ yr⁻¹ or 3.4×10⁴ m³ s⁻¹) to the northern side of the GOM. The remainder comes from other US rivers, Mexico and Cuba (Nipper et al., 2004).

Figure 1 Study sites and sampling areas in the Gulf of Mexico and the Korean Peninsula. Panel (a) shows the sampling area within the northern Gulf of Mexico. Flow in the Mississippi–Atchafalaya river system is split 30 % to the Atchafalaya River and 70 % to the Mississippi River. The box is the sampling area. Panel (b) shows station positions from March 2005. Note that MCH project data are widely distributed across the region. Red, grey and blue stations correspond to sub-regions A (near the Mississippi River), B (between the Mississippi and Atchafalaya) and C (near the Atchafalaya), respectively. Panel (c) shows the sampling area off the west coast of the Korean Peninsula. Panel (d) shows all of the station positions.

At the Old River Control Structure on the lower Mississippi River approximately 25 % of the Mississippi River's water is diverted into the Atchafalaya River, where it mixes with the water in the Red River. The flow in the Atchafalaya River totals 30 % of the total MARS flow (Fig. 1a). Several projects have investigated the relationship between nutrients and the marine ecosystem, and how this leads to hypoxia in the GOM (e.g., Bianchi et al., 2010; Diaz and Rosenberg, 1995, 2008; Forrest et al., 2011; Hetland and DiMarco, 2008; Laurent et al., 2012; Quigg et al., 2011; Rabalais and Smith, 1995; Rabalais et al., 2007; Rabalais and Turner, 2001; Rowe and Chapman, 2002). Strong stratification due to the high freshwater discharge from the MARS, local topography (DiMarco et al., 2010), wind direction and high nitrate concentration all affect hypoxia formation, with upwelling-favorable wind facilitating its development (Feng et al., 2012, 2014).

Figure 2 The Rowe and Chapman three-zone hypothesis, which describes the physical and biochemical processes that initiate and sustain hypoxia on the Louisiana–Texas shelf (Rowe and Chapman, 2002). RMEPs are reduced metabolic end products. Reprinted with permission of Gulf of Mexico Science.

In the northern GOM, the major factor controlling coastal productivity is riverine N input. Rowe and Chapman (2002) defined three theoretical zones over the LATEX shelf close to the Mississippi and Atchafalaya river mouths to predict the effects of nutrient loading on hypoxia along the river plumes and over the shelf. They named these the brown, green and blue zones (Fig. 2). Nearest the river mouths is a “brown” zone, where the nutrient concentrations are high, but the discharge of sediment from the river reduces light penetration and limits primary productivity within the plume. Further away from the river plume is a stratified “green” zone with available light and nutrients that result in high productivity. In this region, the rapid depletion of nutrients is due to biological uptake processes that depend on the season and river flow (Bode and Dortch, 1996; Dortch and Whitledge, 1992; Lohrenz et al., 1999; Turner and Rabalais, 1994). Still further offshore, and also along the river plume to the west, there is the so-called “blue” zone, defined arbitrarily by nitrate concentrations of 1 µM or less, which is dominated by intense seasonal stratification and a strong pycnocline, so that in the surface layer nutrients are limiting at this distance from the rivers and most primary production is fueled by recycled nutrients (Dortch and Whitledge, 1992). It is important to note that RC02 makes clear that the edges of the zones (geographical regimes) are not static, but change over time depending on season, river flow and biological processes (Fig. 2).

The western CSK forms the eastern side of another semi-enclosed basin (the Yellow Sea) and is affected by freshwater discharge from river plumes in the same way as the coastal GOM, although the freshwater flow is considerably less. The Yellow Sea covers about 380 000 km² area with an average water depth of 44 m, and numerous islands are located on its eastern side (Liu et al., 2003). Our specific study area is the mid-western coastal region from the Taean Peninsula to Gomso Bay (Fig. 1c and d).

There is a strong tidal front in the coastal area near the Taean Peninsula due to sea floor topography and the coastal configuration (Park, 2017; Park et al., 2017). The region also contains several bays (Garolim Bay, Gomso Bay and Cheonsu Bay), and is affected by discharges from a large artificial lake (Saemangeum Lake) as well as the freshwater discharge from the Keum River plume that contains high concentrations of nutrients (Lim et al., 2008). Conditions in the mid-western CSK near the Taean Peninsula are similar to the coastal GOM, because of mixing of two different water masses from Gyunggi Bay (Han River) and the Keum River (Choi et al., 1998, 1999). The annual mean flow rates within the Keum River were about 70 m³ s⁻¹ (normal period) and 170 m³ s⁻¹ (flood period; Yang and Ahn, 2008). Precipitation within the catchment was 1208 mm yr⁻¹ during 2003 to 2005 (Yang and Ahn, 2008).

Unlike the coastal GOM, the CSK has increased nitrogen inputs from atmospheric nitrogen deposition (AN-D, which is approximately 5 times higher than in the GOM, Table 2; J. Y. Kim et al., 2010; Luo et al., 2014; Shou et al., 2018; Zhao et al., 2015) and nutrient inputs from the groundwater discharge (J. S. Kim et al., 2010; G. Kim et al., 2011). AN-D has increased in the CSK owing to industrial development in China during the last few decades, which has led to increased atmospheric N emission.

3 Data and methods

3.1 Riverine N data

Hydrographic data from the MCH (Mechanisms Controlling Hypoxia – MCH Atlas) projects in the Gulf of Mexico were collected from the National Oceanographic Data Center (https://www.nodc.noaa.gov, last access: 1 May 2017) covering the period from 2004 through 2007 (Table 1). We excluded cruises MCH M6 and M7 because the threat of hurricanes led to sampling stations in different areas from the other cruises. The study sites and sampling areas are shown in Fig. 1b. Quality control removed inconsistencies and anomalies in the data (e.g., removing outliers, missing data found by linear interpolation). Hydrographic data from the CSK (nutrients, salinity, oxygen) were collected during several cruises (Table 1 and Fig. 1c and d), and the data were put through similar QA/QC routines.

Table 1 Sampling dates for data from Gulf of Mexico projects and the coastal sea of the Korean Peninsula. Winter data are listed for the Gulf of Mexico cruises.

3.2 Atmospheric nitrogen deposition (AN-D) data

AN-D data from around the USA are sparse (Table 2). Most US data have been collected along the east coast of the USA, and the only data in the GOM region were collected near Corpus Christi (∼1 g m⁻² yr⁻¹; Wade and Sweet, 2008). Considerable AN-D could be expected, however, from the large number of petrochemical and fertilizer plants in southern Texas, especially near Houston and along the Mississippi. While there are more data from the Yellow Sea (J. Y. Kim et al., 2010; Luo et al., 2014; Shou et al., 2018; Zhao et al., 2015), they are still limited owing to the broad sampling coverage. While AN-D data in the Asian region were up to 14 g m⁻² yr⁻¹, data from the eastern side of the USA were under 1 g m⁻² yr⁻¹, even lower than in the GOM, suggesting there is currently not a large contribution from AN-D to total N loads in the North Atlantic Ocean. The approximate order of magnitude difference in AN-D concentrations between the GOM and the CSK is due to the continuing industrial development in East Asia and the resulting N emissions (Wang et al., 2016; Zhao et al., 2015). Lamarque et al. (2013) reported model results that cover our study regions, and their model appears to underestimate AN-D at the sampling sites compared with observational data in the GOM (Wade and Sweet, 2008). However, the pattern of AN-D inputs between GOM and CSK from Lamarque et al. (2013) shows around 5 times difference between the two regions, which agrees with our data. Thus, in our model, we used observational data for both regions, as shown in Table 2.

Table 2 Atmospheric nitrogen deposition (AN-D) in the USA and in the Yellow Sea.

3.3 Methodology: N mass balance model

Our model consists of three sub-regions based on sampling locations during MCH cruises (Fig. 3), each of which contains a series of 0.25^∘ square boxes, as followed by Belabbassi (2006). The 0.25^∘ boxes in this study were separated into an upper box and a lower box, based on pycnocline depth, as defined by a sharp change in density and which coincides generally with a minimum change in oxygen concentration of 22.33 µM. We assume steady state conditions, and estimate potential production, which we count as an estimate of potential carbon flux (Fig. 3a). Primary production (PP) above the pycnocline is expected to be higher than below it (Anderson, 1969; Sigman and Hain, 2012), which means that the two layers have different production rates. The difference in PP between upper and lower boxes also depends on the freshwater discharge rate, which determines nutrient input to the upper layer, seasonal variability, and transfer processes between the layers. While chlorophyll can be found below the pycnocline (DiMarco and Zimmerle, 2017), the fact that it is typically associated with low oxygen concentrations suggests that the phytoplankton are either inactive or, more likely, producing at a very slow rate.

Figure 3 (a) Input (blue) and output (red) sources for each 0.25^∘ box in the GOM and CSK (see text for details). (b) Area of each sub-region (red) and boxes affected by direct riverine input (blue) in the GOM. Export N (mixing) represents the advective transport term. The processes of biogeochemical and transport processes of both regions are the same and each in- or output factor is the same in the GOM and CSK. Note that transfer between boxes occurs in both directions alongshore and onshore/offshore and is not a one-dimensional process as suggested in the diagram.

The N mass balance box model is modified from previous models to calculate the net removal of dissolved inorganic nitrogen (DIN) inside each box, which represents potential primary production (PPP; De Boer et al., 2010; G. Kim et al., 2011; Eq. 1). In this model, DIN concentration includes ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$), nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) and nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$).

$$\begin{array}{}\text{(1)}& F_{\mathrm{River}}^{\mathrm{DIN}}+F_{\mathrm{Atmo}}^{\mathrm{DIN}}+F_{\mathrm{Bott}}^{\mathrm{DIN}}-F_{\mathrm{Export}}^{\mathrm{DIN}}-F_{\mathrm{Deni}}^{\mathrm{DIN}}=F_{\mathrm{Removal}}^{\mathrm{DIN}},\end{array}$$

where $F_{\mathrm{River}}^{\mathrm{DIN}}$, an input term, is DIN flux from each river discharge and is calculated with $C_{\mathrm{Box}}^{\mathrm{DIN}}$, the DIN concentration in each box; A_Bott is the bottom area of each 0.25^∘ box; and F_River is the river discharge rate ($C_{\mathrm{Box}}^{\mathrm{DIN}}\times A_{\mathrm{Bott}}\times F_{\mathrm{River}}$). As another input term, $F_{\mathrm{Atmo}}^{\mathrm{DIN}}$ is the flux from atmospheric nitrogen deposition. $F_{\mathrm{Bott}}^{\mathrm{DIN}}$, the benthic flux, is an additional input term in the sub-pycnocline layer box. The 0.25^∘ blue boxes located closest to the Mississippi and Atchafalaya river mouths were assumed to be the only ones affected by riverine input (Fig. 3b). As an output term, $F_{\mathrm{Export}}^{\mathrm{DIN}}$ as an advection term was calculated from the current velocity in each region from observations (Nowlin et al., 1998a, b) and from literature data (Jacob et al., 2000; Lim et al., 2008) and the exchange between boxes from the residence time in each box. Note that water and nutrient exchange can take place through all four sides of each box, so the array is two-dimensional. $F_{\mathrm{Export}}^{\mathrm{DIN}}$ for water mixing was calculated from these factors; $C_{\mathrm{EX}}^{\mathrm{DIN}}$ is the difference in DIN concentration between adjacent boxes; V_S is the water volume of each box; and λ_Mix is the mixing rate of each box $\left(C_{\mathrm{EX}}^{\mathrm{DIN}}\times V_S\times {\mathit{\lambda }}_{\mathrm{Mix}}\right)$. We used a reciprocal of the water residence time that we considered to represent horizontal mixing, i.e., dispersion. Another output term is $F_{\mathrm{Deni}}^{\mathrm{DIN}}$, the denitrification process from the water column, and $F_{\mathrm{Removal}}^{\mathrm{DIN}}$ is removal by biological production. The details of the model definitions are given below in Table 3 and shown in Fig. 3. Each arrow indicates input (blue) and output (red) terms (Fig. 3). Input/output terms vary based on whether the boxes are above/below the pycnocline, while there are separate inputs from the Mississippi and Atchafalaya rivers in the GOM and Keum and Han rivers in the CSK, respectively.

Table 3 Definitions and values used in N mass balance model to calculate DIN removal by biological production.

^a Each 0.25^∘ box. ^b Wade and Sweet (2008) for GOM region. ^c McCarthy et al. (2015). ^d Lee et al. (2012). ^e McCarthy et al. (2015). ^f Qureshi (1995). ^g $F_{\mathrm{Atmo}}^{\mathrm{DIN}}$ of CSK region is used as mean values of Asia data in Table 2, which is initially 5 times higher than that of GOM (1.4×10⁵ mol d⁻¹). ^h The unit of $F_{\mathrm{Sink}}^{\mathrm{DIN}}$ was converted to mol per day from the unit of original data (gN m⁻² d⁻¹) with the area of box (0.25 m × 0.25 m) and molar mass of N (14 g mol⁻¹). All units were converted to mol per day multiplied by the area of the box (0.25 m × 0.25 m).

In order to calculate the net removal of DIN in a box above the pycnocline layer, we used our N mass balance model in Eq. (2).

$$\begin{array}{}\text{(2)}& F_{\mathrm{River}}^{\mathrm{DIN}}+F_{\mathrm{Atmo}}^{\mathrm{DIN}}-F_{\mathrm{Export}}^{\mathrm{DIN}}-F_{\mathrm{Sink}}^{\mathrm{DIN}}=F_{\mathrm{Removal}}^{\mathrm{DIN}}\end{array}$$

The boxes above the pycnocline layer have two input terms: (1) $F_{\mathrm{River}}^{\mathrm{DIN}}$, riverine N, which affects only a subset of boxes along the edge of each region, and (2) $F_{\mathrm{Atmo}}^{\mathrm{DIN}}$, atmospheric nitrogen deposition (AN-D), which affects every box equally. The mean value of Asian data, as shown in Table 2 (J. Y. Kim et al., 2010; Luo et al., 2014; Shou et al., 2018; Zhao et al., 2015), is used for $F_{\mathrm{Atmo}}^{\mathrm{DIN}}$ of the CSK region, which is initially 5 times higher than that of the GOM (1.4×10⁵ mol d⁻¹; Wade and Sweet, 2008). We also considered vertical sinking as an input for the sub-pycnocline layer box and as an output from the upper layer. Other possible input factors might be upwelling/downwelling processes; however, these factors are neglected in the model because both regions are shallow and close to shore (Feng et al., 2014; Lim et al., 2008) and we have no observational data on upwelling/downwelling rates. The output terms are the following: (1) $F_{\mathrm{Export}}^{\mathrm{DIN}}$, the exchange rate between each box (obtained from the different N concentrations in each box and the mass transfer between them), and (2) $F_{\mathrm{Sink}}^{\mathrm{DIN}}$, removal by biological production, including sinking (assuming that any other removal factors are neglected above the pycnocline). We tested the RC02 three-zone hypothesis in the upper box layer, in which we can also examine the horizontal influence (horizontal extent) of the river plume based on production rates.

Below the pycnocline layer we used the revised Eq. (3).

$$\begin{array}{}\text{(3)}& F_{\mathrm{Bott}}^{\mathrm{DIN}}+F_{\mathrm{Sink}}^{\mathrm{DIN}}-F_{\mathrm{Export}}^{\mathrm{DIN}}-F_{\mathrm{Deni}}^{\mathrm{DIN}}=F_{\mathrm{Removal}}^{\mathrm{DIN}}\end{array}$$

Equation (3) has two separate input terms: (1) the benthic flux $F_{\mathrm{Bott}}^{\mathrm{DIN}}$ term contains all the potential input from the bottom sediment (defined here as net DIN release from the bottom sediment) including nutrient regeneration by bacteria, groundwater nutrient inputs and an uptake of nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) and nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$) mainly by sedimentary denitrification (McCarthy et al., 2015; Nunnally et al., 2014), and (2) $F_{\mathrm{Sink}}^{\mathrm{DIN}}$ term as vertical sinking from the box above the pycnocline layer, for which we used data from Qureshi (1995). The unit of $F_{\mathrm{Sink}}^{\mathrm{DIN}}$ was converted to mol d⁻¹ from the unit of original data (gN m⁻² d⁻¹) with area of box (0.25 m × 0.25 m) and molar mass of N (14 g mol⁻¹).

In the GOM, benthic sediments provide excess ammonium to overlying water by regeneration processes such as remineralization (Lehrter et al., 2012; Nunnally et al., 2014; Rowe et al., 2002). Generally, there is an uptake of nitrate and nitrite mainly by sedimentary denitrification (McCarthy et al., 2015) or dissimilatory nitrate reduction to ammonium (DNRA) and assimilation by benthic microalgae (Christensen et al., 2000; Dalsgaard, 2003; Thornton et al., 2007). Due to this, net DIN flux was used as the value of $F_{\mathrm{Bott}}^{\mathrm{DIN}}$, which shows DIN release from bottom sediments to the overlying water column. For example, in the GOM, the sums of nitrate and nitrite fluxes to bottom sediments (e.g., May: −10.05, July: −61.9, August: −48.42 µmol N m⁻² h⁻¹) were similar or smaller than the flux of ammonium from bottom sediments (e.g., May: 203, July: 152, August: 156 µmol N m⁻² h⁻¹) off Terrebonne Bay (McCarthy et al., 2015). In the CSK, the sum of the nitrate and nitrite flux to bottom sediments and ammonium flux are 0.5–1.4 mmol N m⁻² d⁻¹ and 1.3–9.6 mmol N m⁻² d⁻¹, respectively, which indicated that excess ammonium with additional nitrate and nitrite were released from sediments in this region (Lee et al., 2012). The release of nitrate and nitrite in the CSK unlike the GOM can be estimated due to high inputs of nitrogen by groundwater in the CSK (G. Kim et al., 2011) even though there is minor uptake of nitrate and nitrite. Diffusion from groundwater can probably be ignored in the GOM as Rabalais et al. (2002) reported that the groundwater discharge is very low in coastal Louisiana, but is likely important elsewhere and is known to be important in the CSK. Based on this, we averaged and sum the fluxes data of nitrate, nitrite and ammonium from McCarthy et al. (2015) for the GOM and Lee et al. (2012) for the CSK, respectively, and then applied $F_{\mathrm{Bott}}^{\mathrm{DIN}}$ value as 1.2 mmol N m⁻² d⁻¹ in the GOM and 6.2 mmol N m⁻² d⁻¹ in the CSK. Thus, in Eq. (3), the benthic flux term is calculated from existing literature results after considering all DIN fluxes as above (Lee et al., 2012; McCarthy et al., 2015), and then multiplied by the area of each box.

The output terms are (1) $F_{\mathrm{Export}}^{\mathrm{DIN}}$, the exchange rate between each box in the lower layer, and (2) $F_{\mathrm{Deni}}^{\mathrm{DIN}}$, the denitrification rate from the water column. Due to high stratification at the pycnocline, upward transfer of dissolved material from the lower layer to the upper layer is assumed not to occur in our model. Also, denitrification from the water column below the pycnocline is a significant N removal process, which removes up to a maximum 68 % of total N input from the Mississippi River (MR) in the GOM (McCarthy et al., 2015). As the value of $F_{\mathrm{Deni}}^{\mathrm{DIN}}$ in the GOM, we used a direct measurement of denitrification rates from the McCarthy et al. (2015) in the water column (88 µmol m⁻² h⁻¹, which converted to 2.1 mmol N m⁻² d⁻¹) where the stations were exactly same as our sub-regions A, B and C. We assumed this applied only below the pycnocline where oxygen concentrations decrease. However, in the CSK, there is no water column denitrification data because the dissolved oxygen concentration has never been below about 4 mg L⁻¹ during our data periods. Based on this, we estimated that there is very little water column denitrification in the CSK, so we did not count this term in the CSK. Thus, we only considered the sedimentary denitrification term for the CSK region.

Figure 4 Mean ocean current velocities (a) and standard deviations (b) for biweekly periods from August 1993 through December 1994 based on data from LATEX project. Positive values of U show eastward flow; positive values of V show northward flow.

Water transport in the region is generally from the east, i.e., from near the Mississippi River in sub-region A to the west, near the Atchafalaya River in sub-region C during non-summer periods. During summer, the winds change direction from easterly to westerly, blocking the water flow to the west (Cho et al., 1998). We calculated advection from current meter data collected during the LATEX program (Nowlin et al., 1998a, b) from April 1992 to December 1994, from which we determined U (west to east flow) and V (south to north flow) components (cm s⁻¹). Figure 4 shows the mean values of coastal ocean current velocities. The annual range of the currents is 0 to 30 cm s⁻¹ for the longshore component, with standard deviation of about 8 cm s⁻¹, and 0 to 7 cm s⁻¹ for the cross-shelf component, with a similar standard deviation, but these current velocities are not constant and change depending on time and day. The annual current velocities in the CSK are more affected by tidal exchange and the presence of the Yellow Sea Current, but velocities are similar to those in the GOM (Jacob et al., 2000; Lim et al., 2008). The annual range of the currents is around 0 to 28 and 0 to 7 cm s⁻¹ for the cross-shelf component. Thus, we used the mean value of the current velocity for the time of year during each cruise in both the GOM and the CSK for calculating the advective flow in both alongshore and onshore/offshore directions.

To run the box model, we assumed three factors: (1) the study area is in a steady state condition, with equal input sources and outputs, (2) AN-D is evenly distributed across each area and (3) DIN is fully utilized by phytoplankton growth in the layer above the pycnocline, so we can neglect other removal factors. However, in the layer below the pycnocline, as we mentioned above, denitrification, which leads to a main loss of DIN as nitrogen gas, is considered as another output term in Eq. (3). Because we assumed that all DIN removed is fully consumed by primary production above the pycnocline, we can calculate potential carbon fluxes and oxygen consumption using the Redfield ratio (C : N :−O₂ : P = 106 : 16 : 138 : 1). The PPP can be compared with ¹⁴C measurement data (Lohrenz et al., 1998, 1999; Redalje et al., 1994; Quigg et al., 2011) and dissolved oxygen data from MCH mooring C at 29^∘ N, 92^∘ W (4 March to 10 July 2005; Bianchi et al., 2010).

4 Results

4.1 An N mass balance model for the Louisiana–Texas shelf

The existence of the three zones suggested by RC02 has been verified from winter data using nutrient/salinity relationships (Kim, 2018). Figure 5 shows the contour graph based on the mean concentration of DIN at each station during the MCH M4 (March 2005) cruise. For operational and modeling purposes, stations were grouped into three sub-regions – near the Mississippi (A), near the Atchafalaya (C) and an intermediate region (B) between ∼90 and ∼91^∘ W. During summer, it is hard to use nutrient/salinity relationships directly because riverine nutrient inputs are lower and phytoplankton growth causes rapid nutrient consumption over the shelf, leading to low overall nutrient surface concentrations. We calculated the mean DIN in each box, and then used the relationship between DIN and salinity to define the edges of the three zones. Near the coast salinity was consistently low, with high turbidity from the river water discharge. This was labeled the brown (river) zone.

Figure 5 Extent of the three zones defined by RC02 based on the mean concentration of nutrient (DIN) at each station during the MCH M4 cruise in March 2005, showing their correspondence with the three sub-regions used in the box model. Red, grey and blue stations correspond to sub-regions A (near the Mississippi River), B (between the Mississippi and Atchafalaya) and C (near the Atchafalaya), respectively.

A range of N input values from various sources were used in the N mass balance model to estimate PPP and carbon fluxes in the coastal GOM. The PPP rates were highest near the river mouth and we set the boundaries of production for each zone based on our N mass balance model results and mean DIN data. We defined the brown zone as having the PPP rate of over 2 gC m⁻² d⁻¹ because of the high input of N from the river, AN-D and benthic fluxes, and the rate in the blue zone is less than 0.1 gC m⁻² d⁻¹. The PPP rate in the green zone is then between 0.1 and 2 gC m⁻² d⁻¹. Basically, these PPP ranges were set based on synthesized measured ranges of coastal GOM primary production, as defined for near, middle, and far fields of the coastal GOM (Dagg and Breed, 2003; Lohrenz et al., 1999). Note that our model results of the PPP might overestimate the actual production because of light limitation, following RC02.

Figure 6 Areal distributions of the three zones using data from above the pycnocline (a–f), based on N mass balance model results. Colors and numbers represent boxes found in each of the three zones in terms of potential productivity (unit: gC m⁻² d⁻¹). As previously, these distributions use data from below the pycnocline (g–l).

The edges of the three zones above and below the pycnocline layer, based on our N mass balance model results, are shown in Fig. 6a and b. The patterns of the boundaries above and below the pycnocline differ from the edges of the zones. The brown zone was found above the pycnocline on all cruises close to the Mississippi River mouth because of the high nutrient concentrations, but only appeared off the Atchafalaya River in March 2005 (MCH M4). However, below the pycnocline it was found only in April 2004 (MCH M1) in sub-region A. This suggests that vertical transport across the pycnocline rapidly removes the high levels of suspended material that cause light limitation above the pycnocline. In the green zones, the nutrient source is mostly supported directly by the river, with minor additional sources of N from vertical sinking, AN-D and benthic fluxes. We utilized the vertical sinking flux from the sediment trap data from Qureshi (1995) below the pycnocline layer to estimate PPP. This varied between 0.1 and 1.0 gN m⁻² d⁻¹ (Table 3). Typically, in the blue zone where biological production is low, vertical sinking followed by local decomposition is assumed to be the major factor that changes the nutrient concentration in the lower layer. The blue zone is always more extensive below the pycnocline than above it, which suggests there is little or no sub-pycnocline production except close to the coast and/or the river mouths, and reinforces the assumption that any chlorophyll below the pycnocline is inactive (Fig. 6b). Thus, we can identify the horizontal influence of the river plume in the layer below the pycnocline and the variation in the boundaries of the three zones, based on the observed nutrient data from a bottom layer and our N mass balance model. The model suggests that regions of moderate potential productivity extend offshore at least as far as 28^∘30^′ N in sub-region B, both above and below the pycnocline.

4.2 An N mass balance model calibration

The model calibration was done with historic literature data. The literature data suggest that observed PP rates in the green and brown zones of the coastal GOM vary between 0.4 gC m⁻² d⁻¹ (winter) and ∼8 gC m⁻² d⁻¹ (summer; Dagg et al., 2007; Lohrenz et al., 1998, 1999; Redalje et al., 1994). Recently, Quigg et al. (2011) determined the integrated PP rates with ¹⁴C measurements during 2004 in the coastal GOM. The highest integrated PP rates were found near the Mississippi River at 3.5 gC m⁻² d⁻¹ (in July), and near the Atchafalaya River at 2.7–5.9 gC m⁻² d⁻¹ (in May to July; in the brown and green zones). However, the lowest integrated PP rates were on the outer part of the LATEX shelf (the blue zone) at 0.07 gC m⁻² d⁻¹ (in March), 0.04–0.15 gC m⁻² d⁻¹ (in May) and 0.33–0.91 gC m⁻² d⁻¹ (in July). Additionally, Quigg et al. (2011) pointed out that these higher PP values were affected by high riverine nutrients input from the MR that flows westward during that time period.

The actual PP ranges were similar to our model-based PPP (Fig. 6). However, this was different from RC02's brown zone. This might be due to the differences between methods such as ¹⁴C, our N mass balance model and RC02's theoretical model. Typically, RC02 assumed that the brown zone is light-limited due to high sediment turbidity, but our model does not account for this and only considers DIN concentrations. Except for this, our PPP results are similar to direct productivity measurements from the ¹⁴C incubations (Quigg et al., 2011). Our model result (PPP) showed the same range of values as ¹⁴C incubations (e.g., Dagg et al., 2007; Lohrenz et al., 1998, 1999; Quigg et al., 2011; Redalje et al., 1994) in the three sub-regions.

Note that our model assumed all the biological uptake could be converted directly to production rates, which we considered as PPP. The PPP from cruises MCH M1–M8 for samples from above the pycnocline calculated using our model is reasonable based on comparison with previous PP values (Fig. 6a). The PPP ranges (0.01–5.05 gC m⁻² d⁻¹) were similar to previous ¹⁴C measurement PP values of between 0.04 and 5.9 gC m⁻² d⁻¹.

Based on our model calculation, which assumes all the nutrients are available for production, the PPP showed maxima at all times in sub-region A (near the Mississippi River) and minima in sub-region B (between the Mississippi and Atchafalaya rivers), except for MCH M2 in June 2004, when sub-region C had the lowest PPP (Fig. 6a). The high values in sub-region A are due largely to under-utilization of nutrients in regions of high turbidity. As the water flows west under the influence of the Coriolis effect, PPP is expected to decrease as a result of declining nutrient concentrations because of dilution and nutrient uptake during biological production while the water flows to sub-region B. In sub-region C, MCH M4 (March 2005) had the highest PPP among the all MCH cruises. This probably depended on high nutrient concentrations being present during the winter period, when the region was affected by Atchafalaya River nutrient input.

4.3 Model scenarios in the Gulf of Mexico (GOM)

We tested the sensitivity of the model to changes in input/output parameters such as increasing AN-D and decreasing riverine N input. Assuming the model is robust, we investigated three model scenarios based on the nutrient distributions seen during the MCH1 cruise (note that using data from other cruises gives very similar results). In the first scenario, we cut riverine N input 60 % and increased the AN-D input by a factor of 2 based on increasing N emission predictions (Duce et al., 2008; He et al., 2010; Kanakidou et al., 2016; T. Kim et al., 2011; Lawrence et al., 2000; Paerl et al., 2002). In the second scenario, we doubled the amount of AN-D as in scenario 1 and decreased riverine N input by 30 % based on the hypoxia management plan goal (Gulf Hypoxia Action Plan Report, 2001, 2008; Rabalais et al., 2009). In the third scenario, we increased riverine N input by 20 %, assuming the failure of the hypoxia management plan, while we set the AN-D amount equal with the first and second scenarios. Based on our N mass balance model calculation and model scenarios, we can initially estimate carbon fluxes from our PPP rate, and, using the Redfield carbon-to-oxygen stoichiometry ratio (106:138), the overall oxygen balance within the coastal GOM (Table 4).

Table 4 Simulation results for selected model scenarios based on MCH M1 (5–7 April 2004), as described in the text. Biological production is calculated using our N mass balance model, while oxygen demand is calculated by the Redfield stoichiometry ratio (C: $-{\mathrm{O}}_{\mathrm{2}}=\mathrm{106}:\mathrm{138}$; unit: gC m⁻² d⁻¹).

As can be seen in the scenario results for MCH M1 data (Table 4), the riverine N input source is still the major controlling factor in the coastal GOM region even when its contribution is greatly reduced and the AN-D source is doubled. For instance, if we fail to reduce riverine N input in the future (scenario 3), the potential carbon fluxes will increase by 17 % relative to current conditions. In contrast, the AN-D input source only increased to a maximum of 5 % of the total input term; this indicates that AN-D input is still a minor factor in the GOM. If the production is increased, overall oxygen demand will also be increased. The MCH M1 scenario result indicated that the overall oxygen demand would increase approximately 21 % if we fail to reduce riverine N input, likely increasing considerably the area of the hypoxia.

4.4 An N mass balance model in the CSK

As we have done in the GOM, we used our N mass balance model to estimate the PPP in the CSK and define the three different zones (Fig. 7). Similar to the GOM region, the PPP rates were highest near the river mouth, and we set the boundaries of each zone based on our N mass balance model results. Based on nutrient data, as was done for the GOM, we defined the brown zone as having a PPP rate above 1.5 gC m⁻² d⁻¹ because of the increased N sources from the river, AN-D, and the sediment flux. We defined the green zone as having PPP rates between 0.3 and 1.5 gC m⁻² d⁻¹ and the blue zone as having rates of less than 0.3 gC m⁻² d⁻¹.

Figure 7 The distribution of the three zones of the mid-western CSK above the pycnocline (a–d) based on the RC02 hypothesis applied to the N mass balance model. Colors and numbers represent boxes found in each of the three zones in terms of potential productivity (unit: gC m⁻² d⁻¹). As for above, using data from below the pycnocline (e–h).

The seasonal results shown in Fig. 7a and b show that the boundaries of the three zones above and below the pycnocline layer were roughly consistent with the main change coming in summer (August), which is the wet season and sees the highest river discharge. The large size of the green zone in all seasons suggests that AN-D is consistently adding extra nitrogen to the surface ocean along with the riverine N input. This is supported by the fact that the PPP in the blue zone is an order of magnitude higher than for the GOM. Around 90 % of the grid cells in the CSK are in the same zones above and below the pycnocline (Fig. 7a and b) during all four cruises; however, in the GOM (Fig. 6a and b) this was found for fewer than half of the grid cells. This is probably due to the difference in freshwater discharge rate in the two regions, which leads to the much larger stratified area in the GOM than in the CSK.

One question that has not been investigated is the temperature dependence of primary productivity in the two areas. While the GOM is temperate throughout the year, winter temperatures in the CSK fall to ∼5 ^∘C. However, according to the ocean color remote sensing images from near the CSK river mouth reported by Son et al. (2005), primary production in the CSK does not appear to be strongly affected by temperature. The PPP results of our model (0.2 to 2.2 gC m⁻² d⁻¹) agreed with their ocean color remote sensing results (0.4 to 1.6 gC m⁻² d⁻¹) in the CSK. Also, during all seasons, the Keum River consistently supplies high amounts of DIN (average: < 60 µM; Lim et al., 2008) to the coastal zone (especially close to the Keum mouth). We believe, therefore, that the higher value of PPP in winter near the Keum mouth (brown zone in Fig. 7a), is reasonable.

The AN-D input source comes mainly from the Chinese side of the East China Sea (ECS) and this affects the boundaries of the green and blue zones above the pycnocline as it is deposited uniformly across the region. There is also nutrient input from offshore, as the Yellow Sea Bottom Cold Water Mass can up-well during the mixing process and is assumed to supply additional nutrients to the outer shelf (Lim et al., 2008).

4.5 Model scenarios in mid-western CSK

AN-D is currently considerably more important (by approximately an order of magnitude) in the CSK than in the GOM, and it is anticipated that AN-D will likely be a major controlling factor here in the future (Duce et al., 2008; He et al., 2010; T. Kim et al., 2011; Lawrence et al., 2000; Paerl et al., 2002). Because of the lack of research on potential hypoxia scenarios in the Korean Peninsula, we used the same three scenarios in the CSK as were used for the GOM. Similar to GOM results, riverine N input remains the major controlling factor; however, in this area, the AN-D source is more critical than in the GOM region (Table 5). The AN-D input source increased from 20 % to 47 % of the total input under scenario 1, while based on our scenario 3 results, increases in the AN-D input source and riverine N input together will affect biological production by increasing carbon fluxes up to 25 % and oxygen demand up to 32 % if we fail to reduce N input in future (Table 5).

Table 5 Simulation results for selected model scenarios based on CSK (February 2008) data. Biological production is calculated by our N mass balance model. Oxygen demand is calculated by the Redfield stoichiometry ratio (C: $-{\mathrm{O}}_{\mathrm{2}}=\mathrm{106}:\mathrm{138}$; unit: gC m⁻² d⁻¹).

5 Discussion

Most previous model studies in the GOM were focused on predicting the hypoxia area (Bierman et al., 1994; Fennel et al., 2011, 2013; Justic et al., 1996, 2002, 2003; Scavia et al., 2004). For example, Justic et al. (1996, 2003) used a two-layer model incorporating vertical oxygen data, from one station (LUMCON station C6; 28.867^∘ N, 90.483^∘ W), to predict the size of the hypoxia area. Similarly, Fennel et al. (2011, 2013) used their more complex simulation model, which included oxygen concentration as well as a plankton model from Fasham et al. (1990), to predict the size of the hypoxia region in the GOM. Our N mass balance model, in contrast, uses historical data from the LATEX shelf to estimate potential carbon fluxes in the GOM, and calculate the overall oxygen demand from those carbon fluxes. While this affects the total area subject to hypoxia it does not estimate the size of the hypoxic zone.

Figure 8 Distribution of the three zones during cruises MCH M1–M3 based on salinity data (Lahiry, 2007). Areas identified as brown, green or blue zones are shaded accordingly.

In contrast to our model, traditional predictive models have also ignored different nitrogen input sources such as AN-D and SGD. While this is probably reasonable on the Louisiana–Texas shelf, where riverine inputs dominate, it may not apply in other coastal regions. As a result, model studies in this region have concluded that reducing riverine N input is the only solution to decrease the size of the hypoxia area in the GOM (Gulf Hypoxia Action Plan Report, 2001, 2008; Rabalais et al., 2009; Scavia et al., 2013). According to our model results, AN-D is still a minor controlling factor in the GOM; however, in the CSK, the AN-D contributed more to the total nitrogen budget and may be a major controlling factor in the future. This indicates that AN-D should be considered as another input term for nutrient managements, especially in Asia or in other regions where high concentrations are expected. Similarly, nitrogen input from either sediment fluxes or groundwater also need to be considered.

Our zonal boundaries can be compared with the results of Lahiry (2007), who used salinity to define the edges of each zone for the three cruises MCH M1, M2 and M3 (Fig. 8) and defined the edges of the RC02 zones in the coastal GOM based solely on salinity. Lahiry's limited simulation results indicated similar patterns to our model based on DIN concentration near the Mississippi River mouth (e.g., during MCH M1, M2 and M3). Mixing was more conservative in this region than further west because the low salinity water with high nutrient concentrations was less diluted with offshore water.

Away from the MR in sub-regions B and C, however, Lahiry's results gave very different boundaries for the three zones compared with our results (Fig. 8). In particular, the results near the Atchafalaya River were very different (compare Figs. 6 and 8). For example, our data showed only green and blue zones off Atchafalaya Bay during MCH M1, with no brown zone. Similarly, the extent of the blue zones in sub-region C during MCH M2 and M3 is also very different. We believe that our N-model-based classification can cover more complex biological processes than the Lahiry (2007) method, which considers only advection and mixing, and that our N-model is a more sensible way to look at biological processes in the GOM.

Our results also agree with previous studies that demonstrated that both the GOM and CSK regions are N-limited for most of the year (Lim et al., 2008; Turner and Rabalais, 2013). This compares with the results of Sylvan et al. (2007), who reported that the coastal GOM could be P-limited in the MR delta mouth area where our brown zone is located, while RC02 suggested light limitation rather than N- or P-limitation. However, these P-limited conditions appear to occur when N concentrations are very high. In particular, the N∕P ratios in the both the GOM and CSK during our sampling were less than 16, indicating that both regions were N-limited, although a few stations in the brown zone near the MR area had ratios of between 16 and 18 (Fig. 9). These higher N-to-P ratios may result from the high sediment turbidity causing light-limited conditions in this zone near the river mouth (Rowe and Chapman, 2002).

Figure 9 Dissolved inorganic nitrogen (DIN) against dissolved inorganic phosphorus (DIP) during sampling periods in GOM and mid-western CSK. Nearly all samples had an N:P ratio of < 16, which indicated potential N-limited conditions. At a few points near the brown zone the ratio was between 16 and 18; this is where light limitation is expected according to RC02.

It should be remembered, however, that the arithmetic N:P value per se is unimportant in determining nutrient limitation. As long as both nutrients can be measured, it is theoretically possible for phytoplankton to continue to grow. The MARS has generally such an excess of N relative to P that N:P ratios ≫16 can be expected as P concentrations fall, but this does not necessarily mean that productivity is limited, and we never found P concentrations of zero in any of our sub-regions; the lowest P concentration measured during all cruises in the GOM and CSK was 0.2 µM.

Both the GOM and CSK regions receive nitrogen inputs from AN-D, rivers and benthic fluxes. These different nitrogen input sources control coastal productivity, and this may reflect the different nitrogen cycling in the two regions. In the GOM, the riverine N input source consistently dominates coastal productivity and eutrophication, while, in the CSK, AN-D is also becoming a critical controlling factor. In the CSK, however, there is strong tidal mixing of freshwater from the Keum River and/or Gyunggi Bay with nearby coastal water, which results in a tidal front along the offshore region and off the Taean Peninsula during spring and summer. It is this physical mixing that mostly controls the spatial distribution patterns of nutrients and salinity here, particularly below the pycnocline (Lim et al., 2008). The brown zone in the upper layer in the CSK (August 2008) changed to a green zone region below the pycnocline layer as a result of the strong coastal tidal mixing.

RC02 considered their model to be theoretical. In the brown zone, close to the river mouth, they assumed turbidity leads to light-limited conditions. Their results agree well with measured ¹⁴C PP numbers from Quigg et al. (2011), who found the lowest integrated PP is near the MR delta mouth. However, our N mass balance model did not consider light limitation and therefore PPP in the brown zone is high. Such good agreement suggests that our model can be applied to a wide region, while ¹⁴C measurements are typically conducted at a few specific points, as long as such limitations are taken into account.

In the CSK, most previous production studies focused on inshore areas such as estuaries or rivers. Our research focused for the first time on the coastal ocean off the Korean Peninsula. Our results suggest that diverse nitrogen sources need to be recognized as potential issues for future nutrient management concerned with hypoxia, eutrophication and other environmental issues. The agreement between our results and the pattern of production based on satellite-sensing in the CSK (Son et al., 2005) suggests that our model is reasonable.

The results of our changing scenarios represent how the biological processes in these coastal regions may vary as individual nutrient sources change. While our model cannot predict the area of the hypoxic zone, we can investigate the effects of potential flux changes of each factor, such as AN-D, riverine input or benthic fluxes, and calculate the effects of changes in each on PPP and on the overall oxygen balance for the region. We have only considered different input terms of our N mass balance model; output terms such as water mixing rates and the residence time for each box need more detailed study in future work to calculate more realistic production changes in each box.

6 Conclusions

The model suggests that the three-zone theory of RC02 can be applied not only in the northern GOM but also in the CSK region and that three zones can be distinguished based on their nutrient concentration. As a result, we believe that using our N mass balance model to separate different zones based on RC02 may be appropriate not only for large-scale regions like the GOM and CSK but also at small scales such as river or estuary systems. The model also estimates potential primary production and carbon flux based on the inclusion of AN-D data that have not been considered previously (e.g., Bierman et al., 1994; T. Kim et al., 2011). Our results agree well with previous ¹⁴C measurements in the GOM (Quigg et al., 2011) and ocean color remote sensing in the CSK (Son et al., 2005).

Based on CSK cruise data results, we can initially determine where the three different zones are in the CSK. We evaluated our model and tested its sensitivity based on three different scenarios. Through our scenario results, we assume that the AN-D is a considerable factor in the CSK as well as the riverine N input from the Keum River. Reducing nutrient input from the river is critical for hypoxia management policy (Gulf Hypoxia Action Plan Report, 2001, 2008; Rabalais et al., 2009). In addition, these model scenarios will be helpful for future coastal nutrient management and hypoxia management studies in the CSK, especially as AN-D sources become more important.