Dissolved inorganic nitrogen and particulate organic nitrogen budget in the Yucatán shelf: driving mechanisms through a physical–biogeochemical coupled model

Continental shelves are the most productive areas in the seas with the strongest implications for global nitrogen cycling. The Yucatán shelf (YS) is the largest shelf in the Gulf of Mexico (GoM); however, its nitrogen budget has not been quantified. This is largely due to the lack of significant spatio-temporal in situ measurements and the complexity of the shelf dynamics, including coastal upwelling, coastal-trapped waves (CTWs), and influence of the Yucatán Current (YC) via bottom Ekman transport and dynamic uplift. In this paper, we investigate and quantify the nitrogen budget of dissolved inorganic nitrogen (DIN) and particulate organic nitrogen (PON) in the YS using a 9-year output from a coupled physical–biogeochemical model of the GoM. The sum of DIN and PON is here referred to as total nitrogen (TN). Results indicate that the main entrance of DIN is through its southern (continental) and eastern margins. The TN is then advected to the deep oligotrophic Bay of Campeche and central GoM. It is also shown that the inner shelf (bounded by the 50 m isobath) is “efficient” in terms of TN, since all DIN imported into this shelf is consumed by the phytoplankton. Submarine groundwater discharges (SGDs) contribute 20 % of the TN, while denitrification removes up to 53 % of TN that enters into the inner shelf. The high-frequency variability of the TN fluxes in the southern margin is modulated by fluxes from the YC due to enhanced bottom Ekman transport when the YC leans against the shelf break (250 m isobath) on the eastern margin. This current–topography interaction can help to maintain the upwelling of Cape Catoche, uplifting nutrient-rich water into the euphotic layer. The export of TN at both western and northwestern margins is modulated by CTWs with a mean period of about 10 d in agreement with recent observational and modelling studies.

1 Introduction

Continental shelves are the most productive areas in the ocean, widely recognized to play a critical role in the global cycling of nitrogen and carbon (e.g. Fennel, 2010; Liu et al., 2010) with direct implications for human activities, such as fisheries, tourism, and marine resources (Zhang et al., 2019).

The importance of nitrogen budgets in shelves has motivated numerous observational and modelling studies of different shelves in the world (e.g. Fennel et al., 2006; Xue et al., 2013; Ding et al., 2019; Zhang et al., 2019). Their significance lies in that dissolved inorganic nitrogen (DIN) supply fuels primary productivity which in turn impacts the socio-economical and recreational activities in those regions. Furthermore, the exchange of DIN and particulate organic nitrogen (PON) between the shelf and the deep ocean influences the carbon cycle (Huthnance, 1995), and it is strongly correlated with other shelf processes such as acidification, eutrophication, red tides, hypoxia–anoxia zones, pCO₂, and sediment denitrification (Fennel et al., 2006; Seitzinger et al., 2006; Enriquez et al., 2010).

In the GoM (Fig. 1a), with a horizontal extension of almost 250 km, the Yucatán shelf (YS) (Fig. 1b and c) is one of the largest shelves in the world. It has 340 km of littoral extension, representing 3.1 % of Mexico's littoral zone. The Yucatán state in Mexico occupies the 12th place in volume catches and the sixth place in production value of fisheries in the country. The fishery production is increasing every year with a growth of 72 % from 2008 to 2017 (Anuario de Pesca 2017, 2017).

Figure 1 Bathymetry (h_m, m) of the whole model domain. Isobaths: 50, 250, 1000, 2000, 3000, and 4000 m are also shown in grey contours. The small box at the upper right corner of (a) shows the north, east, south, and west boundaries used to compute the inner- and outer-shelf TN cross-shelf fluxes. The yellow dashed box delimits the study area of the Yucatán shelf, where (b) is the surface temporally averaged velocity field (U, m s⁻¹) with magnitude in colour and vectors representing the direction, and (c) is the surface mean kinetic energy (MKE, cm² s⁻²) computed for the year 2010. The smallest yellow box in (a) shows the “notch” area (see text) and the three yellow lines are the mooring locations for transects YUC, PN, and PE. Labels help identify the Deep Gulf of Mexico, Campeche Bay, and Caribbean Sea regions in (a). The inner and outer Yucatán shelf are shown in (c).

Total nitrogen (here, TN is equal to DIN plus PON) fluxes are intrinsically related with the productivity and nitrogen cycling of the shelves. However, sources and sinks of TN are highly uncertain and difficult to quantify. This is partly due to the large spatial and temporal variability associated with the cross-shelf and along-shelf regional nutrient budgets and the difficulty to measure them. Biogeochemical coupled modelling systems are a useful tool to quantify the shelf–open ocean TN exchange, taking into account the different spatial and temporal scales involved in the biogeochemical cycle (Walsh et al., 1989; Fennel et al., 2006; Hermann et al., 2009; Xue et al., 2013; Damien et al., 2018; Zhang et al., 2019).

The physical mechanisms that drive and modulate the cross-shelf transport of nutrients and biogenic material are also poorly known. Shelves are rich dynamical areas in which several processes can coexist at different spatio-temporal scales. Ekman divergence, coastal-trapped waves (CTWs), current interactions with the shelf break, mesoscale structures, vertical mixing, and topographic interactions, among others, are processes that may uplift nutrient-rich waters from the deep ocean into the photic zone of continental shelves (e.g. Cochrane, 1966; Merino, 1997; Roughan and Middleton, 2002, 2004; Hermann et al., 2009; Shaeffer et al., 2014; Jouanno et al., 2018).

In this regard, the YS is a complex system due to the coexistence of different physical processes relevant in its dynamics. One of the first studies in the area is that of Merino (1997), who reported the uplift of nutrient-rich Caribbean waters from 220 to 250 m deep, reaching the YS at the “notch area” (small yellow box in Fig. 1), likely due to the interaction of the Yucatán Current (YC) with the YS. The zonal Caribbean Current of the Cayman Sea turns northwards when reaching the Yucatán Peninsula, forming the strong western boundary YC that flows through the Yucatán Channel, located between the eastern slope of the YS and northwestern Cuba (see yellow line in Fig. 1a). Once inside the GoM, the YC becomes the Loop Current (LC) (Candela et al., 2002), which interacts with the slope of the YS on its eastern side (Cochrane, 1966; Merino, 1997; Ochoa et al., 2001; Sheinbaum et al., 2002), favouring the outcrop of deep nutrient-rich waters to shallower layers over the shelf. However, the mechanisms responsible for this upwelling and its variability remain unclear.

The wind pattern over the YS is characterized by the trade winds (easterly winds) throughout the year, with recurrent northerly wind events during autumn and winter caused by cold atmospheric fronts with relatively short durations (Gutierrez-de Velasco and Winant, 1996; Enriquez et al., 2013). The easterly winds drive a westward circulation over the inner shelf (Enriquez et al., 2010; Ruiz-Castillo et al., 2016). They are also responsible for the upwelling along the zonal Yucatán coast due to divergent Ekman transport (Fig. 2). This upwelling is present year-round along the north and northeast coast of the YS, with intensifications from late spring to autumn (Zavala-Hidalgo et al., 2006).

Figure 2 Seasonal climatology of surface chlorophyll (mg Chl m⁻³) given by the biogeochemical coupled model for (a) winter (January, February, March), (b) spring (April, May, June), (c) summer (July, August, September), and (d) autumn (October, November, December), for the 2002–2010 period. Dashed boxes in (a) denote the three areas in which the validation with observations (black dots) was carried out, i.e. inner shelf, outer shelf, and the upwelling region close to Cape Catoche.

Besides the wind-induced upwelling near the coast, there is also upwelling produced by the interaction of the YC with the eastern YS which is considered the principal mechanism that brings deep nutrient-rich waters over the YS. Observational studies suggest high intrusions of upwelled waters during spring and summer which are suppressed during autumn–winter (Merino, 1997; Enriquez et al., 2013). This seasonal variability is not easy to explain since the YC near the YS does not show such a clear seasonal signal and is dominated by higher-frequency mesoscale variations (Sheinbaum et al., 2016), so several mechanisms have been proposed to understand it. For example, Reyes-Mendoza et al. (2016) show how northerly winds can suppress the upwelling at Cape Catoche. Since these cold front northerly winds are active during autumn–winter, they could explain in part the seasonality of the cold water intrusions. But other mechanisms appear to be important too: CTWs (Jouanno et al., 2016), topographic features and bottom Ekman transport (Cochrane, 1968; Jouanno et al., 2018), extension and intensity of the Loop Current (Sheinbaum et al., 2016), and encroachment and separation of the YC and LC from the shelf (Jouanno et al., 2018; Varela et al., 2018). External (off-shelf) sea level conditions may also generate pressure gradients that oppose the upwelling and explain its seasonality (Zavala-Hidalgo et al., 2006).

Regarding freshwater inflow, a significant source to the YS is related to submarine groundwater discharge (SGD) due to the karstic geological formation of the Yucatán Peninsula (Pope et al., 1991; Gallardo and Marui, 2006), coastal lagoons (Herrera-Silveira et al., 2004), and springs (Valle-Levinson et al., 2011). Due to the complexity of mechanisms and scarcity of observations, the total discharge of SGD into the YS is not well known.

Coupled hydrodynamic–biogeochemical models can be used to establish the TN routes in the marine environment (Fennel et al., 2006). Xue et al. (2013) proposed the first model for TN dynamics in the GoM shelves but excluding the YS. To the best of the authors' knowledge, there are no studies describing the nutrient flux pathways in the YS, so the present work represents the first attempt at a quantitative analysis to understand the biogeochemical cycles and their modulation by physical process in one of the most important socio-economical areas of the southern GoM.

We use a coupled physical–biogeochemical model of the whole GoM to study the nitrogen budget in the YS. The biogeochemical cycles of the YS are poorly known in the GoM and controversies remain regarding its physical dynamics besides the long-term undersampling of biogeochemical variables (Zavala-Hidalgo et al., 2014; Damien et al., 2018), as well as the presence of SGD with unknown fluxes. The main objectives of this study include (i) quantification of the TN budget within the inner and outer YS, (ii) investigation of the sources and sinks of nitrogen in the continental shelf, and (iii) analysis of the physical mechanisms that modulate the cross-shelf TN transport.

2 Model set-up and observational data

2.1 Physical model

The physical model is a GoM configuration of the Regional Ocean Modeling System (ROMS), which is a hydrostatic primitive equations model that uses orthogonal curvilinear coordinates in the horizontal and terrain following (σ) coordinates in the vertical (Haidvogel and Beckmann, 1999). A full description of the model numerics can be found in Shchepetkin and McWilliams (2005) and Shchepetkin and McWilliams (2009). Horizontal grid resolution is ∼5 km, with 36 modified σ layers in the vertical. We used a new vertical stretching option (Azevedo Correia de Souza et al., 2015) that allows higher resolution near the surface. The numerical domain, which covers the whole GoM, is shown in the bathymetry map in Fig. 1a. The model was run for 20 years (1993 to 2012), from which we use 9 years (2002 to 2010) in the present analysis in order to be time-consistent with observational satellite data.

The bathymetry is provided by a combination of the “General Bathymetric Chart of the Oceans” (GEBCO) database (http://www.gebco.net/, last access: April 2015) with data collected during several cruises in the GoM. The initial and open boundary conditions for temperature, salinity, and velocity come from the GLORYS2V3 reanalysis which contains daily averaged fields (Ferry et al., 2012). The model is also forced with hourly tides obtained from the Oregon State University TOPEX/Poseidon Global Inverse Solution (TPXO) (Egbert and Erofeeva, 2002). Hourly atmospheric forcing comes from the “Climate Forecast System Reanalysis” (CFSR) (Dee et al., 2014). These include cloud cover, 10 m winds, sea level atmospheric pressure, incident short- and long-wave radiation, latent and sensible heat fluxes, and air temperature and humidity at 2 m. These variables are provided at ≈38 km horizontal resolution and are used to estimate surface heat fluxes in the model using bulk formulae (Fairall et al., 2003). The model uses a recursive three-dimensional MPDATA advection scheme for tracers, a third-order upwind advection scheme for momentum, and a turbulence closure scheme for vertical mixing from Mellor and Yamada (1982).

2.2 Biogeochemical model

The biogeochemical model is described in Fennel et al. (2006) and is based on the Fasham et al. (1990) model which takes nitrogen-based nutrients as a limiting factor. The model is solved for seven state variables, namely nitrate (NO₃), ammonium (NH₄), phytoplankton (Phy), zooplankton (Zoo), chlorophyll (Chl), and two pools of detritus: large detritus (LDet) and small detritus (SDet). Details of the model algorithm and coupling to ROMS can be found in Fennel et al. (2006). An important aspect of this model is a better simulation of denitrification processes at the sediment–ocean interface in the bottom of the continental shelves.

Initial and boundary conditions for the biogeochemical variables were obtained from an annual climatology of NO₃, NH₄, and Chl. The climatology was calculated using all available profiles with the highest quality control from the World Ocean Database (Boyer et al., 2013) and profiles obtained from the XIXIMI cruises carried out by CICESE. The DIVA optimal interpolation (Troupin et al., 2012) scheme was used to interpolate the individual profiles in the climatology to the model grid. DIVA takes into account the coastline geometry, sub-basins, and advection to reduce errors due to artefacts in the interpolation.

The XIXIMI cruises provided profiles of nutrients and chlorophyll in the southern GoM, which helps to reduce the bias between the northern and southern parts of the GoM. The cruises encompass the region between 12 and 26^∘ N and −85 and −97^∘ W and were carried out within the scope of the “Consorcio de Investigación del Golfo de México” (CIGoM) project (“Gulf of Mexico Research Consortium” project in English).

Close inspection of the shelf dynamics through maps of the temporally averaged velocity field $U=\stackrel{\mathrm{‾}}{u}$, $\stackrel{\mathrm{‾}}{v})$ (Fig. 1b), where the overline denotes the temporal mean, and mean kinetic energy MKE = $\mathrm{0.5}({\stackrel{\mathrm{‾}}{u}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v}}^{\mathrm{2}})$ (Fig. 1c) allows us to delimit the shelf into two areas. The first is the inner shelf, delimited by the 50 m isobath where the strongest YS velocities develop (Fig. 1b) and where most of the MKE is enclosed (Fig. 1c). The second area is the outer shelf between the 50 m and the 250 m isobaths, with the latter isobath representing the shelf break.

The TN examined in this study is taken as the sum of DIN and PON, with DIN = NO₃ + NH₄, and PON = Phy + Zoo + SDet + LDet (Xue et al., 2013). The cross-shelf nitrogen fluxes are calculated as

$$\begin{array}{}\text{(1)}& Q_{\mathrm{50}\mathrm{m},\mathrm{250}\mathrm{m}}=\underset{-\mathrm{50},-\mathrm{250}}{\overset{\mathit{\eta }}{\int }}{u}_{\mathrm{cross}}\left(N\right(\mathrm{d}z,\end{array}$$

where u_cross is the velocity component normal to the 50 or 250 m isobaths, η is the model sea level, and N can be any component of the TN. The TN cross-shelf fluxes are computed for the north, east, south, and west boundaries for both the inner and outer shelf, as indicated in the inset in the upper right corner of Fig. 1a. Accordingly, the total budget is obtained as the integral over the area of the shelf and over the depth of the water column for both the inner and outer shelves. The budget also includes the loss to denitrification and to burial in the sediments, which are taken into account for the quantification of the TN budget as sinks of nitrogen.

The initial concentration and boundary conditions at the edges of the GoM model domain (Fig. 1a) of the biogeochemical variables NH₄, Phy, Zoo, Chl, and pools of detritus are set to a small and positive value of 0.1 mmol N m⁻³ following Fennel et al. (2006, 2011) and Xue et al. (2013). As mentioned in these references, the model quickly adjusts internally to proper variable values within days to weeks. Moreover, these boundaries are far away from the YS and therefore the fluxes across the inner and outer YS determined internally in the model are not impacted by possible inconsistencies at the GoM open boundaries. Given the lack of data for Mexican rivers and groundwater fluxes, the same approach is followed for freshwater inputs as also followed by Xue et al. (2013). The biological model parameters used in this study are those shown in Table 1 of Fennel et al. (2006), except for the vertical sinking rates which were reduced about 10 %, to fit the depth of the deep chlorophyll maximum (DCM) observed with the APEX profiling floats (see Fig. A4). The model does not include an explicit compartment for nitrogen in the form of DON, although it can be included as in the work of Druon et al. (2010), which adds semi-labile DOC and DON as state variables to the original Fennel et al. (2006) model. They comment on the difficulties of validating the model with observations and highlight open questions even in the definition of both DOC and DON pools (see also Anderson et al., 2015). Considering these difficulties and uncertainties, our approach is to use, initially, more basic models to understand their capabilities and build and employ more comprehensive ones later on; so the inclusion of DON and/or DOC compartments is left for future studies.

2.3 Freshwater sources

Two riverine systems account for 80 % of the freshwater discharge into the GoM, the Mississippi–Atchafalaya system with 18,000 m³ s⁻¹ and the Usumacinta–Grijalva system with 4500 m³ s⁻¹ (Dunn, 1996; Yáñez Arancibia and Day, 2004; Kemp et al., 2016) (see Appendix B). Freshwater contributions to water volume, salinity, temperature, and DIN concentration are included as grid-cell sources into the model. Apart from the two main systems, a total of 81 freshwater sources are included, taking into account freshwater discharges in the Florida, Texas, and Yucatán shelves from the years 1978 to 2015. For the US rivers the daily data were obtained from the US Geological Survey (USGS) (https://www.usgs.gov/, last access: January 2017) and the Gulf of Mexico Coastal Ocean Observing System (GCOOS) (https://products.gcoos.org/, last access: January 2017).

Although the YS has no rivers, freshwater inputs play a key role impacting the local ecosystem (Herrera-Silveira et al., 2002). These inputs come from SGD linked to the “cenotes” ring (sink holes) system inland. The freshwater flux, temperature, salinity, and nutrient concentrations for these sources are not well known. Monthly climatological values were calculated for the Mexican rivers and SGD systems, using temporally scattered information found in the literature (e.g. Rojas-Galaviz et al., 1992; Milliman and Syvitski, 1992; Poot-Delgado et al., 2015; Conan et al., 2016) and a data collection effort within Mexican institutions led by Jorge Zavala-Hidalgo (personal communication, January 2015) and from the GOMEX IV cruise of CINVESTAV (Centro de Investigación y Estudios Avanzados in Merida Yucatan) within the CIGoM project. During this cruise, a total of 71 profiles of NO₃, potential temperature, salinity, and chlorophyll were collected at standard depths from 2 to 20 November 2015. The localization of the profiles is shown in Fig. 2. Therefore, fluxes from US rivers forcing the model present inter-annual variability but Mexican freshwater sources only include a climatology due to lack of information (see Appendix B for more details).

The nitrogen concentration for freshwater sources is essentially DIN. For most of the northern rivers (e.g. Mississippi and Atchafalaya), PON is also considered where available (Fennel et al., 2011; Xue et al., 2013). For the remaining freshwater sources, including the SGD system of YS, the PON contribution is set as a constant small value of 1.5 mmol N m⁻³ due to lack of data.

3 Model evaluation for the YS

The model dynamics and its biogeochemistry are validated to guarantee the simulation is able to reproduce basic features of the observations in the GoM, particularly in the YS. Model statistics including biases of physical and biological variables are computed to have some idea of their impact on the estimation of the TN budget over this shelf. Since this is a basin-scale coupled model, a general evaluation of the results and their statistics is carried out considering sea surface temperature, mixed-layer depth, mean kinetic energy, surface chlorophyll, and deep chlorophyll maximum over the whole Gulf of Mexico with emphasis on the YS. The results are presented in Appendix A.

Figure 3 Seasonal climatology of bottom temperature (^∘C) for (a) spring and (b) autumn, for the period between 2002 and 2010. The corresponding vertical sections, indicated by the zonal black line in (a), are shown for (c) spring and (d) autumn. The contours in (a) and (b) denote the 50 and 250 m isobaths. The black contour in (c) and (d) shows the upwelling isotherm of 22.5 ^∘C.

3.1 YS in situ data comparison

Recall that upwelling into the YS is more intense during spring–summer and weaker in autumn–winter (Ruiz-Castillo et al., 2016; Merino, 1997). While the model presents upwelling during all the simulated months, this seasonal behaviour is represented in the model climatologies shown in Fig. 2. The figure also shows the position of oceanographic stations occupied during the GOMEX IV oceanographic cruise and delimitation of three areas of particular interest: the inner shelf, the outer shelf, and the upwelling region at Cape Catoche. The climatology of the YS bottom temperature (Fig. 3) shows that cold waters enter into the shelf during spring in agreement with the enhancement of chlorophyll concentrations (Fig. 2b). The zonal vertical cross sections show that the isotherm of 22.5 ^∘C, which traces the upwelled water (Cochrane, 1968; Merino, 1997), outcrops into the shelf during spring (Fig. 3c). This is not the case in autumn (Fig. 3d), and the upwelling is weaker (Fig. 2d).

A point-by-point comparison between the model results and the in situ observations is shown using only data for November from 2002 to 2010 in the model, for compatibility with the observation dates (Figs. 4–7). Since the simulation is for different years, we only expect to reproduce basic features of these observations. The range of temperatures at different depths shown by the model agrees well with those observed during GOMEX IV (Fig. 4). The mean temperature of the observations is 25.5±2.9 ^∘C, while the model mean temperature is 24.3±3.7 ^∘C. The bias of −1.3 ^∘C is deemed acceptable considering the model mean is a 9-year mean whereas the mean from observations is from just one month and a different year. A critical area to be evaluated is the upwelling region (see dashed box in Fig. 2a), the bias there is −1.1 ^∘C with a root-mean-square error of 1.68 ^∘C. This means that the model tends to be slightly colder than the observations even inside upwelling waters.

Figure 4 Comparison between in situ data and simulated temperatures (^∘C). Temperature values correspond to each hydrographic station, averaged over three depths; (a) between the surface and 25 m depth, (b) between 25 and 50 m depth, and (c) between 55 and the deepest measured concentration ($z\sim -\mathrm{150}$ m). Black dots correspond to the observed values and open grey circles to the simulation. Vertical grey lines are the temporal standard deviation of the simulated values, as these are temporally averaged over all Novembers from 2002 to 2010. Vertical black lines delimit the group of stations for the inner shelf, outer shelf, and upwelling area.

The model mean salinity is 36.5±0.2, which matches the 36.5±0.2 from observations (Fig. 5). Whilst surface salinity in the model is in relatively good agreement with observations (Fig. 5a), differences become more important at deeper layers (Fig. 5b and c). The root-mean-square error of model salinity (0.23) as well as the bias (−0.04) is low, which tends to underestimate the salinity observations. These low differences are also found in the bias for the upwelling area, although the model overestimates the salinity by 0.21 there. The model is able to represent the main characteristics of the Caribbean Subtropical Underwater coming from the Caribbean Sea (Merino, 1997) and the Gulf Common Water from the GoM (e.g. Enriquez et al., 2013) within the YS. The warm and high-salinity Yucatán Sea Water at the surface described in Enriquez et al. (2013) is present in the model too, although temperatures do not exceed 31^∘ (not shown) as in observations.

Figure 5 Same as Fig. 4, but for salinity.

For Chl, the model results fall within the range obtained from fluorometer observations in the inner shelf, outer shelf, and upwelling areas (Fig. 6). The mean observed Chl (0.52±0.58 mg Chl m⁻³) is slightly larger than the model results (0.44±0.42 mg Chl m⁻³) but within the 1 standard deviation range, with a bias of −0.08 mg Chl m⁻³ and a root-mean-square error of 1.16 mg Chl m⁻³. Notice that there is agreement in Chl concentration between the model and observations in the three layers between 150 m depth and the surface (Fig. 6a–c). In the upwelling area the model has lower concentrations than observations with a bias of −0.39 mg Chl m⁻³ and a root-mean-square error of 1.39 mg Chl m⁻³; although the bias is relatively low, it needs to be taken into consideration for the TN budget. Additionally, a comparison with observed mean chlorophyll vertical profiles over the YS is presented in Fig. A9 of Appendix A. Profiles have a similar structure but the model tends to underestimate the DCM.

Figure 6 Same as Fig. 4, but for chlorophyll concentrations (mg Chl m⁻³).

To evaluate the temporal behaviour of the model Chl, time series of the surface chlorophyll averaged over the shelf are compared to similar time series from satellite surface chlorophyll from MODIS (see Fig. 8c and Appendix A for a description of the satellite product) during the simulation period. Mean values of satellite surface Chl are 0.38±0.09 and 0.36±0.13 mg Chl m⁻³ in the model. Besides reproducing temporal mean and variability of the surface chlorophyll, the model is able to reproduce a positive trend present in the 9 years of satellite data. No trend is present in any of the biogeochemical forcings of the model, and determining which physical mechanisms produce it requires further investigation (see below).

The simulated nutrient concentration depicts values with a similar order of magnitude ($\sim \mathrm{3.1}\pm \mathrm{4.6}$ mmol N m⁻³) as the observed profiles ($\sim \mathrm{3.7}\pm \mathrm{5.2}$ mmol N m⁻³) (Fig. 7). Surface nutrient concentrations are underestimated by a 1.7 (mmol N m⁻³) compared to observed profiles (Fig. 7a). At subsurface depths (25–55 m), the model tends to underestimate the NO₃ concentrations; however, in the upwelling area, model NO₃ concentrations are closer to the observed values with a bias of −0.7 mmol N m⁻³ and larger standard deviations for both the model (4.0±5.0 mmol N m⁻³) and observations (4.81±6.33 mmol N m⁻³) (Fig. 7b). The temporal variability of the modelled NO₃ is larger than the observed NO₃ at the surface and bottom as shown by the largest standard deviation in Fig. 7b. Below 55 m the modelled and observed NO₃ are in good agreement in both the outer shelf and the upwelling area (Fig. 7c). Again, these model results are deemed consistent with observations and are in the range of other values reported in the literature (Merino, 1997). Comparison of similar budgets from other shelves in the GoM can be made (e.g. Xue et al., 2013), though clear interpretation of similarities and differences between them may be difficult given the differences in dynamics and nitrate sources and sinks controlling the budgets on each shelf. One could easily compute budgets per unit area or length for a more sensible comparison among different shelves but in the literature only total budgets are available (see Table 2 in Xue et al., 2013). In that regard, the inner and outer shelf areas are ∼74 and ∼91 km², respectively. The TN concentrations for each shelf can be extracted by averaging over the 9-year simulation the integrated values of Fig. 8, with 1$\mathrm{4.6}\pm \mathrm{0.83}\times {\mathrm{10}}^{\mathrm{16}}$ and $\mathrm{7.42}\pm \mathrm{0.89}\times {\mathrm{10}}^{\mathrm{16}}$ mmol N, for the inner and outer shelf, respectively.

Figure 7 Same as Fig. 4, but for nitrate concentrations (mmol N m⁻³).

Figure 8 Temporal series of TN (thick black line), DIN (thin black line), and DIN (thin grey line) in millimoles of nitrogen, spatially integrated over (a) the inner shelf and (b) the outer shelf. Panel (c) shows the temporal series from monthly satellite chlorophyll (black, mg Chl m⁻³) and from the model outputs (grey) averaged over the whole Yucatán shelf. Dashed thick lines are the trend indicated by the linear fit for the TN or chlorophyll time series, where thinner dashed lines are the respective 95 % confidence intervals. Equations of each linear fit are TN (inner shelf) = 2.33×10¹² d + 4.2×10¹⁶, TN (outer shelf) = 2.40×10¹² d + 7.0×10¹⁶, Chl (satellite) = 0.0010 months + 0.28, and Chl (model) = 0.0010 months + 0.30. Notice that the trend is positive for all the temporal series.

In addition, the model sea level elevation and surface ocean currents are compared against altimeter products in Appendix A, where YC variability and transport from the model are compared with data from three moorings located on the slope close to the eastern YS rim described in Sheinbaum et al. (2016).

4 TN budget and cross-shelf transports in the YS

Time series of spatially averaged TN over the YS suggest a positive trend over the 9 simulated years. The trend is seen in both the inner and the outer shelves (Fig. 8a and b). This, perhaps, could be expected given the positive trend in both model and satellite surface Chl mentioned before (Fig. 8c). Varela et al. (2018) report a cooling trend of the inner YS and suggest it may be associated with an eastward shift of the YC. We searched for possible connections between the trends in chlorophyll and TN and indices measuring the position and strength of the YC in the model but found no correlation. This is an interesting problem currently under investigation and to be reported elsewhere.

In the inner shelf there are similar total integrated values of DIN and PON (Fig. 8a). This indicates the presence of a very “efficient” biogeochemical cycle in the inner shelf (see explanation below). By contrast, in the outer shelf, DIN values are larger than PON (Fig. 8b) probably because the integration in the outer shelf includes a large volume below the euphotic zone. Temporal series of integrated TN show a combination of low-frequency variability associated with the seasonal cycle as well as interannual variability, but longer period integrations are required to properly investigate the latter.

To understand the high TN variability in the YS, quantification of the cross-shelf fluxes becomes necessary. Their impact on the TN budget and the physical mechanisms modulating such fluxes are investigated next.

Cross-shelf fluxes are quantified for the two compartments, the inner and outer shelves (Fig. 1b), and for all the boundaries of each compartment. A schematic view of the main incoming and outgoing pathways of cross-shelf TN fluxes is shown in Fig. 9. The yearly averages of the spatially integrated cross-shelf fluxes are shown in Table 1.

Figure 9 Scheme of the TN budget for the Yucatán shelf. Black and white arrows denote cross-shelf flux for the outer and inner shelf, respectively. In blue is the PON, in red the DIN, in green the freshwater DIN sources (Rivers), and in yellow the sinks of TN due to denitrification (DNF). The values are expressed in mol N yr${}^{-\mathrm{1}}\times {\mathrm{10}}^{\mathrm{10}}$. Negative values indicate sink, whereas positive indicates a source of TN. The isobaths that delimit the inner (50 m depth) and outer (250 m depth) shelves are also highlighted.

Table 1 Nutrient budget in moles of nitrogen per year for the inner (50 m isobath) and outer (250 m isobath) Yucatán shelf, computed at each boundary (N, W, E, and S; see Fig. 1a) using cross-shelf velocities. The flux of nutrients is integrated through the water column and temporally averaged using the period 2002–2010 to compute the budget from daily model fields. Positive (negative) values represent sources (sinks) of nutrients. Denitrification is always a nitrogen removal process.

^a Freshwater sources are considered only for the inner shelf. Inner can be taken as a source or sink of nitrogen only for the outer shelf. ^b The positive trend of total nitrogen observed in the temporal series during 9 years is also taken into consideration to close the budget.

For the inner shelf, both PON and DIN are imported through its northern and eastern boundaries and exported through the west and south borders. The inner shelf acts as a source of PON for the Campeche Bay at the southwest margin. The major source of TN for the inner shelf is from the outer shelf via the Cape Catoche upwelling, representing 80 % of the total, while freshwater sources contribute the other 20 %. Although the latter is a relatively large source of nitrogen, its relevance seems to be confined to the NW part of the inner shelf. In general, there is compensation between the DIN and PON concentrations in the inner shelf (Fig. 8a) due to an efficient biogeochemical cycle whereby almost all the DIN imported into the shelf is consumed by the phytoplankton and thus converted into PON. The efficiency relies on the shallowness of the inner shelf (∼50 m depth), because, if strong mixing conditions are present, organic matter will distribute throughout the shallow water column. This is enhanced during winter, when vertical wind-driven mixing and convective processes are strong enough to reach the sea bottom. Additionally, vertical shear likely generated by bottom friction can lead to instabilities and vertical mixing able to break the stratification and carry nutrients to the euphotic zone. During summer months, vertical mixing is weaker (not shown). Turns out that, in the model, vertical velocities in the inner shelf are quite intense and upward throughout the year (∼5 m d⁻¹), carrying nutrients to the euphotic layer. The cause of these vertical velocities is under investigation using a higher-resolution model configuration.

By contrast, in the outer shelf, the largest inputs of PON and DIN are advected from its southeastern corner. The eastern boundary is a source of DIN but a sink of PON for the outer shelf. Therefore, the budget reveals that the PON exported to the inner shelf is produced in the outer shelf and not advected from Caribbean waters. The contribution of TN from the inner shelf to the outer shelf represents only 1.5 % of the total inputs.

Over the outer shelf the fluxes of nutrients and organic matter are driven by a westward wind-driven circulation (Ruiz-Castillo et al., 2016) and exported to the deep GoM and the Campeche Bay through the north and west borders, respectively. This represents a source of DIN, Phy, and Zoo to these oligotrophic regions.

Figure 10 (a) Map of surface chlorophyll (mg Chl m⁻³), averaged over the 9 simulated years. The three characteristic isobaths are denoted. A total of 9 years of averaged cross-shelf fluxes along the 250 m isobath at the western boundary of (b) TN, (c) DIN, and (d) PON (mmol N m⁻² s⁻¹). Negative values indicate westward flux, i.e. TN flux from the shelf to the Campeche Bay. The area delimited by dashed lines shows the location of the filament depicted in (a), in the NW of the YS.

In that regard, the model reveals a quasi-permanent thin filament of Chl that is advected from the northwest corner of the outer shelf to the west of the Campeche Bay (Fig. 10a). A vertical section of the cross-shelf fluxes along the 250 m isobath in the western YS (TN, Fig. 10b) shows that while the export of organic matter to the open sea is concentrated in the surface layers (Fig. 10d), the bottom layer presents a net DIN export (Fig. 10c). The climatological average over 9 years of simulated Chl shows that this filament is intensified during winter (not shown), although it is present during the whole simulation period. Sanvicente-Añorve et al. (2014) studied the larval dispersal for coral reef ecosystems in the southern GoM. They show that the northwestern corner of the outer YS acts as a sink region for larvae. Similar to other coral reef systems, they attributed the sink to the influence of circulation patterns that lead to a unidirectional dispersion pattern during the whole year. Our model results seem to support this idea.

Denitrification is a form of anaerobic microbial respiration in which nitrate and nitrite are finally reduced to molecular nitrogen (N²). It represents a major sink for bioavailable nitrogen. The spatio-temporal average rate of denitrification for the YS is 1.11±0.13 mmol N m⁻² d⁻¹. Our results suggest that denitrification removes up to the 53 % of the TN in the inner shelf, a significant percentage that agrees with estimates from other shelves in the GoM (Xue et al., 2013). On the other hand, denitrification in the outer shelf only removes 9 % of the TN. Our results also indicate that denitrification rates tend to increase with time for both inner and outer shelves (not shown), in concert with TN concentrations (Fig. 8a and b). This is expected since denitrification is a reduction process; hence an increase in nitrate concentration means more available DIN to be reduced to N².

4.1 Physical modulation of cross-shelf TN flux by CTWs

Many physical process coexist at different spatio-temporal scales in the YS that modulate the cross-shelf transport of nutrients and organic matter. We suggest that at least two processes are responsible for such modulation: CTWs and interaction of the YC with the eastern shelf break.

CTWs can be generated by wind forcing over irregular bottom topography along the coast and have been the subject of investigation for a long time (e.g. Clarke, 1977). In the GoM, CTWs are forced by alongshore winds and then travel anticlockwise with the coast on its right until they reach the western portion of the Yucatán Peninsula, mainly associated with cold fronts (Dubranna et al., 2011; Jouanno et al., 2016). CTWs have a signature in sea level that is well captured in relatively high-resolution models such as the one used in the present study (∼5 km). In their modelling study, Jouanno et al. (2018) suggest that CTWs may influence the Yucatán upwelling pulses. In this study, the presence of CTWs is corroborated and its effect on the modulation of cross-shelf nutrient fluxes at the west margin of the YS is exposed.

The presence of CTWs in the model simulations is evidenced in the Hovmöller diagram along the 50 m isobath shown in Fig. 11. Phase speeds are in the range of [2–4] m s⁻¹ in agreement with observations (Dubranna et al., 2011) and other models (Jouanno et al., 2016).

Figure 11 (a) Snapshot of sea level anomaly (η, m) for the simulated year 2005. (b) Hovmöller diagram of η along the 50 m isobath from January to April of the year 2005. Red dots in (a) denote the latitude and longitude shown at the bottom of (b), from Florida to the Yucatán Peninsula.

The cross-shelf TN in the YS western inner-shelf boundary exhibits high-frequency variability. The daily climatology of the wavelet power spectrum of wind stress, sea level anomaly (SLA), and western boundary cross-shelf TN temporal series for the inner shelf of Fig. 12 show that both along-shelf wind stress and changes in SLA may be linked with the cross-shelf TN variability in the inner shelf. The three variables show maximum energy during winter times when CTWs are expected to be more intense, and the wind increases its magnitude due to the incursion of the “Nortes” (cold front winds).

Figure 12 Wavelet power spectrum for time series averaged over the western 50 m isobath of (a) cross-shelf total nitrogen flux vertically integrated (TN, mmol N m⁻¹ s⁻¹), (b)  along-shelf wind stress (τ_along, N m⁻²), and (c) sea level anomaly (SLA, m). The temporal series are detrended, normalized, and filtered by a Lanczos high-pass filter with a cut-off of 15 d.

It is worth mentioning that the wavelet power spectrum for the whole 2002–2010 period (not shown) depicts an interesting intensification of cross-shelf flow (and nutrient fluxes) during 2003, 2004, 2009, and 2010 which coincides with El Niño Modoki events (Ashok et al., 2007; Ashok and Yamagata, 2009). The possibility of such a connection deserves further investigation.

Figure 13 (a) Cross-correlation spectral analysis of time series over the western 50 m isobath, indicating the square coherence coefficient between along-shelf wind stress (τ_along, N m⁻²) and cross-shelf total nitrogen flux vertically integrated (TN_cross, mmol N m⁻¹ s⁻¹) in black and between sea level anomaly (SLA, m) and TN_cross in grey. The black horizontal line indicates the 95 % confidence interval. Analysis for the 9 simulated years based on daily outputs with a 30 d window. Before analysis, the temporal series are detrended, normalized, and filtered by a Lanczos high-pass filter with a cut-off of 15 d. Panel (b) shows the phase or anti-phase in degrees of both coherence analysis of (a) τ_along and TN_cross in black and SLA and TN_cross in grey.

To further examine the relationship between these physical and biogeochemical variables, results from a cross-correlation spectral analysis are shown in Fig. 13 for the time series used in the wavelet analysis of Fig. 12. The variability of along-shelf wind stress and cross-shelf TN fluxes shows significant coherence in the 8–10 d period band at nearly zero phase lag (Fig. 13). Coherence between cross-shelf TN fluxes and SLA is also coherent in the same band (peaks at 8 and 8.4 d) but 180^∘ out of phase. This is consistent with offshore Ekman transport produced by along-shelf northerly winds triggering nutrient and organic matter fluxes across the western boundary of the YS and negative SLA at the coast.

Propagation of CTWs is evident in the Hovmöller diagram of Fig. 11 and most certainly modulates the cross-shelf TN transport. The coherent 8–10 d period band (and also at other higher frequencies, e.g. 5–6 d period) is in agreement with those reported in the literature for CTWs in the GoM (e.g. Jouanno et al., 2016). Since the coherence analysis is carried out here using time series of spatially averaged quantities (from 20^∘30^′ N to almost 22^∘ N, approximately 100 km), possible phase lags are probably masked.

4.2 Influence of the Yucatán Current in the coastal upwelling

Observational studies suggest that favourable-upwelling winds at the northern Yucatán coast are present all year round (Ruiz-Castillo et al., 2016; Pérez-Santos et al., 2010). Cold SSTs on the YS vary seasonally and are particularly characterized by a cold water band on the inner YS very close to the coast that appears in spring and continues until the beginning of autumn (Ruiz-Castillo et al., 2016; Zavala-Hidalgo et al., 2006). Pérez-Santos et al. (2010), using 10 years of sea surface wind data from QuikSCAT, show that Ekman transport is the main contributor to the upwelling over the north YS (93 %), with Ekman pumping only contributing 7 %.

This upwelling regime requires a supply of cold and deeper nutrient-rich waters from the open ocean to maintain the observed biological productivity on the YS.

The main import of TN to the YS is through the southeast and eastern YS boundaries via mechanisms related to the dynamics of the Yucatán and Loop currents and their interaction with the YS shelf break, such as intensification, separation and/or encroachment from the coast, bottom boundary layer transport, advection, instabilities, eddies, and the presence of particular topographic features (e.g. submarine canyons. The reader is referred to Roughan and Middleton (2002) for a discussion of upwelling mechanisms on the East Australian Current that appear to be relevant here too as several local studies indicate (Cochrane, 1966; Merino, 1997; Zavala-Hidalgo et al., 2006; Enriquez et al., 2010, 2013; Enriquez and Mariño Tapia, 2014; Carrillo et al., 2016; Jouanno et al., 2016, 2018).

On the other hand, the export of TN to the deep GoM through the YS northern margin can also be related to advection by the YC and associated mesoscale structures (Roughan and Middleton, 2002; Carrillo et al., 2016; Enriquez and Mariño Tapia, 2014).

Correlation analysis between the strength of the cross-shelf flow from the YC, PON, and DIN fluxes at the eastern margin, all vertically integrated, shows high values at seasonal timescales (Fig. 14). The time series are filtered by a 30 d moving-average window to remove high-frequency variability. The square of the correlation coefficients (r²) for DIN and PON against the vertically integrated YC is indicated at the top of Fig. 14. These results indicate that TN fluxes are well correlated with the strength of the current.

Figure 14 Temporal series for the 9 simulated years of cross-shelf Yucatán Current component (YC), total nitrogen (TN), dissolved inorganic nitrogen (DIN), and particulate organic nitrogen (PON), vertically integrated and averaged over the isobath 250 m of the eastern boundary. The square of the correlation coefficients (r²), between YC and the biogeochemical variables, is shown at the top. The temporal series are filtered by a moving average of 30 d to remove daily variability.

To investigate the possible role of the position and trajectory of the YC and its closeness to the YS in the upwelling (Enriquez et al., 2010, 2013; Jouanno et al., 2018), we computed an index measuring the closeness of the YC core to Cape Catoche and the notch areas in the model, which are two places where water tends to upwell (Merino, 1997; Jouanno et al., 2018). The index depicts no seasonality that could be directly connected to strong(weak) upwelling during spring (autumn). This is an indication that seasonality of the inflow of nutrient-rich water into the YS is probably more influenced by other processes as discussed in Reyes-Mendoza et al. (2016), such as cold front winds that can stop the upwelling or other non-local perturbations.

One of the important mechanisms suggested since Cochrane (1966) to be responsible for the YS eastern boundary upwelling and the nutrient flux towards the coast is bottom Ekman layer transport produced by interaction of the YC with the upper slope and shelf break. The stress exerted by the intense along-shore velocity of the YC on the topography generates an Ekman spiral at the bottom boundary layer and a net depth integrated transport to the left, i.e. a cross-shelf transport towards the shelf in the boundary layer. For example, Shaeffer et al. (2014) using glider observations find that bottom Ekman transport can explain up to the 71 % of the bottom cross-shelf transport variability on the southeastern Australian shelf produced by the East Australian Current.

Figure 15 Flux of TN (Q_TN) computed by the bottom Ekman transport (U_bE, m² s⁻¹) for the 9 simulated years (blue) compared with the bottom-most layer TN flux (grey, mmol N m⁻¹ s⁻¹) over the Ekman bottom layer for (a) temporal averages and (b) spatial averages over the 250 m isobath, where the superimposed black line is the bottom Ekman transport filtered with a 90 d moving average. (c) Vertically integrated TN flux along the eastern 250 m isobath, averaged over latitude and over the 9 simulated years in mmol N m⁻¹ s⁻¹. Shaded areas denote the standard deviation of the averages.

Here we present modelling evidence that bottom Ekman layer transport could be the precursor of the upwelling in Cape Catoche. The bottom Ekman transport (U_bE, m² s⁻¹) can be taken as $U_{\mathrm{bE}}=-{\mathit{\tau }}_{\mathrm{by}}/\left({\mathit{\rho }}_{\mathrm{o}}f\right)$, where τ_by is the bottom stress computed by ${\mathit{\tau }}_{\mathrm{by}}={\mathit{\rho }}_{\mathrm{o}}Cd{v}_{\mathrm{b}}\sqrt{u_{\mathrm{b}}^{\mathrm{2}}+v_{\mathrm{b}}^{\mathrm{2}}}$, with Cd=1×10⁻³ as the drag coefficient, u_b and v_b the bottom velocities at the 250 m isobath, f the Coriolis frequency, and ρ_o=1025 kg m⁻³ the reference potential density of the sea water. The analysis shows that the time-mean U_bE is toward the shelf (defined positive here), and is well correlated with the bottom cross-shelf water transport (r²=0.71, ci = 95 %) calculated directly (Fig. 15a). The Ekman transport is calculated from the theoretical formula (i.e. stress divided by the coriolis frequency), whereas the direct transport is calculated using the bottom velocities and integrating on the last grid cell. We should mention here that the bottom grid cell at this depth has a vertical size of ∼20 m. Using standard formulas to estimate the width of the Ekman layer (e.g. Cushman-Roisin and Beckers, 2011; Perlin et al., 2007) from bottom velocities or stresses and stratification, we obtain values ∼10–30 m. Therefore the layer is not really resolved by the model grid. The correlation is also large over time (r²=0.78, ci = 95 %) as shown in Fig. 15b. Figure 15c shows that the vertically integrated TN transport averaged over 9 simulated years and over latitude is towards the coast at 250 m depth, that is, at the bottom-most model layer which is considered here to be the bottom Ekman layer.

Comparison between bottom-layer Ekman transport and the time mean vertically integrated TN transport across the eastern 250 m isobath indicates that bottom Ekman transport is responsible for 65 % of the TN that is entering the shelf. The mechanisms that explain the remaining flux need to be further investigated and are probably related to meanders, eddies, topographic features, and other processes. Moreover, bottom Ekman transport can be arrested by stratification and may not be dominant everywhere along the YS east coast as has been documented in other western boundary upwelling regions (e.g. Roughan and Middleton, 2002, 2004). Our goal here was only to estimate the size of the TN fluxes related to the bottom Ekman layer and determine its relative importance.

5 Model uncertainties

The bias of the model with respect to observations described in Sect. 3 (see also Appendix A) is analysed in order to examine how uncertainties may impact the TN budget calculation.

The model tends to overestimate (underestimate) Chl (SST) in winter and underestimate (overestimate) Chl (SS) in summer. This bias produces more intense upwelling at Cape Catoche during spring than in summer. In fact, upwelling waters are still present during winter (Fig. 2) but not in observations. The filtered seasonal time series of bottom Ekman transport shown in Fig. 15b (black line) depict the same pattern: they indicate that water from the Caribbean Sea entering into the YS (via bottom Ekman transport) increases during spring towards the summer, decreases during autumn, and increases again during winter. This is in agreement with Fig. 3.

In the water column, the model underestimates NO₃ concentration a maximum of 15 % and is also about 5 % colder than the observed vertical profiles. These biases can impact the growth of phytoplankton whose maximum growth rate (Eppley, 1972) depends on temperature and nutrient concentration (Fennel et al., 2006). However, since phytoplankton only represent 15 % of the TN, the overall impact on TN of these biases is estimated to be less than 3 %. The main point we are trying to make is that, although there are model biases, the main processes that control the TN budget in the YS are well captured by the model simulation, particularly the Cape Catoche upwelling, which together with the southeastern boundary represent the main entrance of TN to the YS.

6 Summary and concluding remarks

The TN budget, the main nutrient transport pathways and their modulation by physical processes over the Yucatán shelf have been investigated using a 9-year simulation from a ROMS physical–biogeochemical coupled model for the GoM. Our work provides a first general view of the shelf physical–biogeochemical coupled system, schematized in Fig. 16.

Figure 16 Schematic view of the main physical processes that modulate the cross-shelf transport of TN in the Yucatán shelf.

The results indicate that TN, especially DIN, enters the Yucatán outer shelf through the southeastern and eastern margins. The TN is then driven by a westward shelf current and then exported to the deep GoM and Campeche Bay through the northern and western boundaries, respectively. In the inner shelf, the biogeochemical nitrogen-based cycle seems to be very efficient for NO₃ remineralization and consumption by the phytoplankton, converting most of the DIN to PON. The freshwater sources represent an important contribution of about 20 % to the DIN concentration, although it is restricted to the northwest of the Yucatán Peninsula. Denitrification represents the main sink of nutrients for the inner shelf, removing more than the 50 % of the nitrogen. The inner shelf contributes to the TN of the outer shelf at its western edge, but this contribution is less than 2 %, indicating weak fluxes from the inner to the outer Yucatán shelf. By contrast, the outer shelf is the main nitrogen supplier of the inner shelf, particularly of PON from the eastern margin. A quasi-constant filament in the outer shelf western border represents an important source of both organic and inorganic nitrogen for the oligotrophic Campeche Bay.

Surface Ekman layer dynamics at the western and northwestern shelf borders play an important role in the transport of nutrient and organic matter to the Campeche Bay and deep central GoM. Part of the high-frequency variability of the TN fluxes at the western YS boundary are correlated and in phase with the along-shelf wind-stress modulating the variability of TN across the western shelf of the Yucatán in the 5–10 period band. These high-frequency TN fluxes are also correlated with changes in SLA at similar periods, which are also typical of CTWs found in the GoM. Coherence is 180^∘ out of phase and consistent with negative SLA resulting from offshore Ekman transport. These exchanges are enhanced during winter due to cold frontal atmospheric systems Nortes.

The advection by the YC dominates the nutrient concentration import to the YS through the southeastern border. This advection, together with the influence of mesoscale structures, control the export of nutrients to the deep GoM at the northern margin. A different process modulates the flux of nutrients at the eastern YS margin. The YC flowing parallel to the slope plays an important role in the intrusion of DIN into the shelf. Initial estimates carried out in this study suggest that, in the model, bottom Ekman layer transport explains the deep TN flux through the eastern YS boundary. There is a positive mean transport (into the shelf) over the 9 simulated years along the eastern shelf break so friction generated between the YC and the shelf break produces a net bottom Ekman transport towards the shelf. This produces the upwelling at Cape Catoche on the eastern shelf, but it is not the only process at work: external and remote forcings appear to control its seasonality (e.g. winds, CTWs); in addition, other upwelling mechanisms such as divergence and convergence from current separation and encroachment and eddy-current interactions with topographic features (e.g. submarine canyons) may be important too and must be analysed in future research.

Appendix A: Model evaluation

A1 Basin-scale model evaluation

This study is focused on the Yucatán shelf region, whose hydrodynamics and biogeochemical outputs were previously validated before the analysis. This appendix provides a summary of this validation to provide evidence that the basic features of the whole GoM circulation are correctly represented in the model. The time-mean eddy kinetic energy (EKE, m² s⁻²) map computed from AVISO geostrophic velocities (http://www.marine.copernicus.eu, last access: June 2018) is used to compare it with the EKE from the model (Fig. A1a and b) for 17 years (1995–2012). The model is able to capture the main features of the eddy field exhibited by the altimeter product as well as the main structure of the Loop Current. In particular, the mean and standard deviation of the eddy kinetic energy field are reasonably captured by the model. The model produces a hook-like pattern of EKE in the western part of the GoM, between 24 and 28^∘ N, that is more evident in the standard deviation of model EKE (Fig. A1c and d). This pattern is not so clear in the AVISO maps but has been identified using Lagrangian data (e.g. Gough et al., 2019). It is associated with the GoM western boundary current that isolates the western continental shelf from the open ocean. EKE is higher in the model, particularly at the Yucatán channel and Florida Straits, probably due to the higher resolution of the model (∼5 km) compared to the altimeter product (∼28 km).

Seasonal climatology of the sea surface temperature (SST) is also compared with the Aqua MODIS satellite products (http://modis-atmos.gsfc.nasa.gov, last access: June 2018) (Fig. A2). The model SST shows a good agreement with seasonal cycles exhibited by the satellite data. The overall bias for the deep GoM SST (depths > 1000 m) is in the range [−0.21, +0.21] ^∘C. Larger differences are found near the coast. The model tends to underestimate the coastal SST during winter and overestimates it during summer. Nevertheless, these differences are less than 0.5 ^∘C, and on average differences are on the order of 0.05 ^∘C with a standard deviation of 0.4 ^∘C. The relatively good agreement between model and data is perhaps not very surprising considering that observed air temperatures are provided to the model to compute heat fluxes using bulk formulae. At the same time, no flux correction is applied in the model so it is important to confirm that there is no drift in the simulation.

In order to evaluate the mixed-layer depth, a total of 2629 ARGO float profiles, available in the period between 1995 and 2012, are compared with the mixed-layer depth given by the model in the deep GoM. This is an important quantity in terms of biogeochemical behaviour since the gulf is an oligotrophic region in which the vertical advection of nutrients controls primary production to the photic layer (Fennel et al., 2006; Xue et al., 2013; Damien et al., 2018). The biogeochemical cycles are partly controlled by the difference between the deep and dark nutrient-rich waters and the upper ocean layer where the availability of light promotes the growth of phytoplankton and hence zooplankton. Figure A3 shows that the model can reproduce the seasonal cycle of the mixed layer in the GoM, with deepening during winter and shallowing during summer seasons (Damien et al., 2018; Portela et al., 2018). The model depicts shallower mixed-layer depths during summer and deeper depths during winter than the Argo observations. The higher variability of the observed data is likely related to mesoscale structures and submesoscale process which can locally deepen the mixed layer or make it more shallow (e.g. Boccaletti et al., 2007; Fox-Kemper et al., 2008; Levy et al., 2012) not fully represented by the model. Despite the differences found, the bias between observations and model mixed-layer depths are on the order of 1.4 m.

The lack of spatio-temporal biological datasets to validate biogeochemical models in the GoM is a well-known problem (Walsh et al., 1989; Damien et al., 2018). Only satellite-derived surface chlorophyll concentration is available with enough spatial and temporal cover but only at the surface. These observations give us a general overview of the chlorophyll temporal and spatial distribution patterns at basin scales. Monthly mean time series (2002–2010) of chlorophyll a concentration from Aqua MODIS and SeaWiFS 9 and 4 km (when available) satellite products are used for a basin-scale model evaluation. The temporal series averaged for the whole deep GoM (i.e. excluding highly productive coastal areas with less than −1000 m depth) show a good agreement between the coupled model and the observations. The model tends to overestimate the Chl in winter and underestimate it in summer (Fig. A5). Despite some exceptional years (e.g. 1999), the modelled chlorophyll concentration values fall in the range exhibited by the satellite products. Mean satellite Chl is 0.1448±0.04 mg Chl m⁻³ in contrast with mean modelling Chl values of 0.1433±0.09 mg Chl m⁻³.

Observations of the vertical chlorophyll structure are available from eight APEX profiling floats with 537 profiles of Chl from 0 to 2000 m every 10 s within the GoM (Pasqueron de Fommervault et al., 2017) (Fig. A4). A more detailed description of this database is provided by Hamilton et al. (2017), and the Chl data calibration is explained in Pasqueron de Fommervault et al. (2017). The resulting profiles give valuable information to evaluate biogeochemical models through the water column, in contrast to the surface-only information from satellite measurements (Pasqueron de Fommervault et al., 2017; Damien et al., 2018). The comparison shows that the model is able to reproduce the depth of the DCM measured by the floats. The DCM seasonal cycle is also well represented by the model. It is interesting to note the high dispersion in the data, revealing the large Chl variability found in the deep GoM.

A2 Regional chlorophyll model evaluation

In addition to the comparison of the surface chlorophyll temporal series with satellite products (Fig. 8c), in situ spatially averaged vertical profiles of chlorophyll from three GOMEX cruises carried out during November 2015, August 2016, and July 2018 are also compared with the model chlorophyll profiles averaged for July, August, and November from 2002 to 2010. The observed profiles superimposed in blue are shown in Fig. A6. The result shows that the model is also able to reproduce the large variability of the observed data. The highest values of chlorophyll from model profiles are found at the surface layers, between 5 and 15 m depth. Values higher than 6 mg Chl m⁻³ represent only 0.64 % of the total simulated points, while for observations, the percentage is about 0.06 %, and these are also located at the surface between 10 and 35 m depth.

A3 Regional altimetry and ocean current comparison

The variance of the absolute dynamic topography (ADT) from AVISO, which is the sea surface height above the geoid obtained as the sum of the sea level anomaly (SLA) and the mean dynamic topography, is compared with the variance of the sea level of the model output (Fig. A7). Observed and model sea surface height (SSH) variance have good resemblance. There are slight differences in the northern coast of the YS. Remember, however, that the accuracy of the altimeter observations is reduced in shallow areas (Vignudelli et al., 2011).

The variability and magnitude of the current over the shelf is also compared against the GlobCurrent product (http://www.globcurrent.org/, last access: January 2019) (Rio et al., 2014). Since the current velocity over the YS is a westward wind-driven flow (Ruiz-Castillo et al., 2016), a comparison with only the geostrophic velocity contribution might not represent the whole state of the velocity field. In this regard, the GlobCurrent product is the result of combining geostrophic altimeter velocity with the addition of the wind-forced Ekman velocity contribution under ocean mixed-layer model assumptions. The results are shown in Fig. A8. The model correctly represents the mean surface current magnitude and direction over the shelf, and the highest differences are found close to the Yucatán coast (Fig. A8a and b). The variability ellipses (Fig. A8c) show that the current variability over 9 years from the model agree with those from observations. Near the northern Yucatán coast, values are lower in both model and data. However, the model ellipses are zonally oriented in contrast to the meridional orientation of the ellipses from the satellite product. The other important difference is found at the west coast of the YS, where the model exhibit southwestward-oriented ellipses whereas the satellite shows a westward orientation. This might influence the direction of the TN fluxes at the west YS boundary, a subject which is addressed in the section Model uncertainties. Similarly to the previous comparison with the AVISO product, significant differences are found near the coast but there are probably significant errors in the data (Vignudelli et al., 2011). In contrast, the YC is well represented by the model in terms of its spatio-temporal variability, although its magnitude is overestimated, which again is probably an effect of better model spatial resolution.

A4 Yucatán Current evaluation

The CICESE-CANEK mooring sections monitoring the flow in the region during 2009–2011 are shown in Fig. 1a (yellow zonal line). The current velocity normal to the three mooring transects shown in Fig. 8 of Sheinbaum et al. (2016) during the years 2009 to 2011 is used for validation. They observe that the YC (YUC transect) was located between the surface and 800 m depth, which agrees with our model results shown in Fig. A9a. The core of the YC is located over the west Yucatán slope, and its mean of 1.18 m s⁻¹ is in very good agreement with observations (Sheinbaum et al., 2002, 2016). The model also shows that the highest standard deviation is at the surface on the western side of the channel with a value of 0.3 m s⁻¹ (Fig. A9d), in contrast with the 0.4 m s⁻¹ found by Sheinbaum et al. (2016). They argue that this variability is due to changes in the current position and the counterflow. At deeper layers (below 900 m), the model shows that the current flows towards the GoM at the centre of the section. On both the western and eastern side of the section, the model is able to reproduce the southward flow as shown in Sheinbaum et al. (2016).

For sections PE and PN (Fig. A9b and c), the model exhibits mean velocities of 0.24±0.24 and 0.36±0.29 m s⁻¹; these values are lower than 0.6±0.7 and 0.4±0.6 m s⁻¹ reported by Sheinbaum et al. (2016) (Fig. A9e and f). One should consider that the model has high variability. Moreover, these sections may or may not be influenced by the core of the Loop Current. Sheinbaum et al. (2016) estimate a reduction of about 30 %–50 % in the maxima of the mean when the Loop Current core is not passing over the section moorings. The southward flow over the slope of section PE below 1000 m is well represented by the model (Fig. A9b), as well as the flow across the whole PN section (Fig. A9c).

Figure A1 Comparison of 17-year (1995–2012) averaged eddy kinetic energy (EKE, m² s⁻²) calculated based on (a) AVISO SSH product and (b) ROMS model simulated SSH. Panels (c) and (d) show the standard deviation for altimeter and model EKE, respectively.

Figure A2 Seasonal climatologies of SST (^∘C) for the GoM (2005 to 2012). Comparison between (a–d) satellite SST product and (e–h) model SST.

Figure A3 (a) Location of the 2629 ARGO profiles used to compute the mixed-layer depth (h_ρ, m). (b) Climatology comparison of mixed-layer depths for ARGO profile floats (black boxes) and the model (grey boxes). Vertical lines in the boxes denote standard deviation.

Figure A4 Seasonal comparison of chlorophyll profiles in milligrams of chlorophyll per cubic metre, taking all the available Apex floats (Pasqueron de Fommervault et al., 2015), in order to evaluate the deep chlorophyll maximum (DCM). Grey dots are the data observed from Apex floats; the average profile is shown in grey. In black is the averaged profile of the model data with its respective standard deviation in dashed black lines.

Figure A5 (a) Temporal series of surface chlorophyll in milligrams of chlorophyll per cubic metre from the satellite and model for the whole deep GoM. Standard deviations from the spatial averages are shown in shaded blue areas for the satellite and dashed black lines for the model. The monthly climatology of the temporal series is shown in the upper part of the figure, where vertical bars indicate standard deviation from the temporal mean. In (b) are represented the correlation coefficients of both monthly temporal series and their respective linear fit as black lines. The slope of the linear fit is 0.25.

Figure A6 Profiles of chlorophyll (mg Chl m⁻³). In black are the model profiles temporally averaged for July, August, and November of the 9 simulated years. Superimposed in blue are the observed profiles of the three GOMEX cruises carried out during November 2015, August 2016, and July 2018.

Figure A7 Variance (σ², m) for the range of the years 2002–2010 of (a) absolute dynamic topography (ADT) extracted from the AVISO altimeter product and (b) mean sea level (η) from the ROMS outputs.

Figure A8 Mean current magnitude and velocity vectors averaged for the years 2002 to 2010 in metres per second. (a) Mean geostrophic plus Ekman currents from the GlobCurrent product. (b) Mean total current from the model; (c) velocity variability ellipses of the GlobCurrent product (black) and model outputs (blue).

Figure A9 Mean model velocity (m s⁻¹) component normal to the three sections, (a) YUC, (b) PE, and (c) PN, depicted in Fig. 1a, for the years 2009 to 2011, and to be compared with Sheinbaum et al. (2016). Positive velocities are in grey (contours every 0.1 m s⁻¹) and negative velocities in white (contours every 0.03 m s⁻¹); the dashed black contour indicates zero velocity. Panels (d)–(f) show the standard deviations for each transect (contours every 0.05 m s⁻¹ for d and e, and every 0.01 m s⁻¹ for f).

Appendix B: Freshwater source inputs

As already described in Sect. 2.3, freshwater, and nitrogen input from 81 major rivers and freshwater systems are included in the coupled model simulation. The Mississippi and Atchafalaya riverine systems are the largest fluvial source in the GoM (red and blue points in Fig. B1a). Their nitrogen delivers tripled in the last decades and are meaningfully correlated with the coastal DIN concentration in the northern GoM (Xue et al., 2013). Nutrient and transport of this system generally peaked in spring-summer in agreement with the time series inputs shown in Fig. B1b.

The Usumacinta–Grijalva rivers system (green points in Fig. B1a) is the most important freshwater source in the southern GoM (Xue et al., 2013). The highest riverine discharge of this system, accompanied by an enhancement of nutrient concentration, occurs during winter and decreases during spring (Fig. B1c). In contrast to the northern riverine sources, the data available for the southern freshwater sources of the GoM are scarce or undersampled. In order to obtain southern freshwater inputs, a time series is built based on the composite of temperature, salinity, volume transport and nutrient concentration at the location of the freshwater sources or near it. The information is obtained from values reported in the literature (Milliman and Syvitski, 1992; Herrera-Silveira et al., 2002, 2004; Yáñez Arancibia and Day, 2004; Hudson et al., 2005; Herrera-Silveira and Morales-Ojeda, 2010; Poot-Delgado et al., 2015; Kemp et al., 2016; Conan et al., 2016) and from observational hydrographic stations near the Yucatán coast (Sect. 2.3). In general, the information reported does not cover a large time series. Therefore, all the information is used to build a climatology which will serve as freshwater model input. An example of this climatological inputs is depicted in Fig. B1c and d for the Usumacinta–Grijalva rivers and the Yucatán freshwater sources (magenta point in Fig. B1a) (Yáñez Arancibia and Day, 2004; Poot-Delgado et al., 2015; Conan et al., 2016). The YS receive freshwater and nutrient inputs from spring and runoff from mangrove areas, lagoons, and cenotes. High nutrient concentrations are reported for YS lagoons (e.g. Dzilam Lagoon) (Herrera-Silveira and Morales-Ojeda, 2010).

As one can notice, the inter-annual variability is visible in northern riverine systems, while is absent in the southern freshwater sources due the lack of information. Moreover, it is essential to note that the small-scale variability in most of the GoM rivers structure is not fully resolved by the horizontal resolution of our model configuration.

Figure B1 (a) Model bathymetry with the location of the input freshwater sources (white points). Shown as red and blue points is the Mississippi and Atchafalaya riverine system, in green is the Usumacinta and Grijalva riverine system, and in violet is the freshwater system of the YS. The panels below show the temporal series of water transport (m³ s) and the DIN (NO₃ + NH₄) fluxes (mmol N s⁻¹) for the systems: (b) Mississippi–Atchafalaya, (c) Usumacinta–Grijalva, and (d) Yucatán shelf freshwater.