Budget of the total nitrogen in the Yucatan Shelf: driving mechanisms through a physical-biogeochemical coupled model

Continental shelves are the most productive areas in the seas with strongest implications for global Total Nitrogen (TN) cycling. The Yucatan shelf is the largest shelf in the Gulf of Mexico (GoM), however, its general TN 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 the Yucatan Current, coastal upwelling, Coastal Trapped Waves (CTWs) and bottom Ekman transport. Through a nine years output of a coupled physical-biogeochemical model of the GoM, the TN budget in the Yucatan shelf 5 is quantified. Results indicate that the main entrance of inorganic nitrogen is through its southern and eastern margins. The TN is then advected to the oligotrophic deep GoM and to the deep Campeche bay. The analysis also shows that the inner shelf (50 m isobath) is efficient in terms of TN, since all the DIN imported into the shelf is consumed by the phytoplankton. Rivers contribute 20 % of the TN, while denitrification removes up to 53 % of TN that enters into the inner shelf. The highfrequency variability of the TN fluxes are modulated by the Yucatan Current in the south and by bottom Ekman transport 10 produced by this current against the shelf-break (250 m isobath) in the east. 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 10 days in agreement with recent observational and modelling studies.

Ekman transport (Enriquez et al., 2010;Reyes-Mendoza et al., 2016). Dynamical processes that may also be important for the upwelling include the extension and the intensity of the Loop Current (Bunge et al., 2002;Jouanno et al., 2018), CTWs (Jouanno et al., 2016), and the bottom Ekman transport. External (off-shelf) sea level may also generate pressure gradients that oppose the upwelling and explain its seasonality.
Regarding the freshwater inflow, a significant amount of sources have been identified in the YS related to groundwater 5 inflows due to the karstic geological formation of Yucatan peninsula (Pope et al., 1191;Gallardo and Marui, 2006), coastal lagoons (Herrera-Silveira et al., 2004), and sinkholes. Due to the complexity of processes and scarcity of observations, the total discharge of freshwater into the YS is not well known.
Owing to the spatial and temporal scarcity of biogeochemical observations, the budget and pathways of nitrogen to the GoM shelves has been poorly quantified. This is specially truth for the southern part of GoM. Coupled models demonstrated to be 10 an efficient tool to establish the routes of the TN (Fennel et al., 2006). In this sense Xue et al. (2013) proposed a first model for the TN dynamics for the GoM shelves, excluding however the YS. From our knowledge, there are no studies describing the nutrient flux pathways in the widest Yucatan shelf. Therefore, the present work represents the first quantitative analysis aiming to understand the biogeochemical cycles and their modulation by physical process at one of the most important socioeconomical areas of the southern Gulf of Mexico. The physical model is the Regional Ocean Modeling System (ROMS) which is a hydrostatic primitive equations model that uses orthogonal curvilinear coordinates in the horizontal and terrain following (sigma) coordinates in the vertical (Haidvogel and Beckmann, 1999). A full description of the model numerics can be found in Shchepetkin and McWilliams (2005) and 20 Shchepetkin and McWilliams (2009). Horizontal grid resolution is ∼ 5 km, with 36 modified sigma 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 Figure 1a. The model was ran for 20 years (1993 to 2012), from which we used 9 years (2002 to 2010) in the present analysis in order to be consistent with the observational data. 25 The bathymetry is provided by a combination of the "General Bathymetric Chart of the Oceans" (GEBCO) database (http://www.gebco.net/) with data collected during several cruises in the GoM. The initial and boundary conditions for temperature, salinity and velocity come from the GLORYS2V3 reanalysis (Ferry et al., 2012). 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 longwave radiation, latent and sensible heat fluxes, and air temperature and humidity 30 at 2 m. These variables are provided at ≈38 km horizontal resolution.

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 limiting factor. The model is solved for seven state variables, namely: Nitrate (NO 3 ), Amonium (NH 4 ), Phytoplankton (P hy), 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). 5 Initial and boundary conditions for the biogeochemical variables were obtained from an annual climatology of NO 3 , NH 4 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 executed by the CICESE group. The DIVA optimal interpolation (Troupin et al., 2012) scheme was used to combine the individual profiles in the climatology using the model grid. DIVA takes into account the coast line geometry, sub-basins and advection to reduce errors due to artifacts in the 10 interpolation.
However, it is well known that the available data density in the GoM is skewed towards the north. The XIXIMI cruises reduce this bias providing profiles of nutrients and chlorophyll in the southern GoM. The cruises encompass the region between 12 • N and 26 • N and -85 • W and -97 • W, and were executed within the scope of the "Consorcio de Investigación del Golfo de México" where most of the MKE is enclosed (Figure 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. 20 The TN examined in this study is taken as the sum of the Dissolved Inorganic Nitrogen (DIN) and the Particulate Organic Nitrogen (PON), with DIN = NO 3 + NH 4 , and PON = P hy + Zoo + SDet + LDet. The cross-shelf nitrogen fluxes are calculated as: where u cross is the velocity component normal to the 50 m or 250 m isobaths, η is the model sea level anomaly and N z can 25 be any component of the TN. Accordingly, the total budget is obtained as: where x, y and σ indicate longitude, latitude and sigma (terrain-following) layer, with n the number of horizontal grid points or vertical model layers.
The initial concentration of the biogeochemical variables (NH 4 , P hy, Zoo, Chl, and pools of detritus) is set to a small and positive value following Fennel et al. (2006Fennel et al. ( , 2011Xue 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 Figure 3).

Freshwater sources 5
Two riverine systems account for 80% of the freshwater discharge into the GoM, the Mississippi/Atchafalaya system with 18,000 m 3 s −1 , and Usumacinta/Grijalva system with 4500 m 3 s −1 (Dunn, 1996;Yáñez Arancibia and Day, 2004;Kemp et al., 2016 -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 by Dr. Jorge Zavala-Hidalgo (personal communication). Therefore the US rivers present inter-annual variability but it is absent in the Mexican rivers.

15
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 underwater groundwater sources linked to the "cenotes" (caves) system. The freshwater flux, temperature, salinity, and nutrient concentrations for these sources are not usually known. For modeling purposes, these are inversely estimated by adjusting the model results to the few temperature and salinity in situ data available and to chlorophyll satellite images. 20 The nitrogen concentration for freshwater sources is essentially DIN. For most of the northern rivers (e.g., Mississippi and Atchafalaya), PON is also considered (Fennel et al., 2011;Xue et al., 2013). For the remaining freshwater sources the PON contribution is set as a constant small value of 0.1 mmolN m −3 .
We also use hydrographic and biogeochemical observational stations taken during the GOMEX IV cruise in the study area.
During this cruise a total of 71 profiles of NO 3 , potential temperature, salinity and chlorophyll were collected at standard The lack of spatio-temporal biological data sets to validate biogeochemical models in the GoM is a well known problem (Walsh et al., 1989;Damien et al., 2018). The most abundant dataset is the satellite derived surface chlorophyll concentration.

30
These observations give us a general overview of the chlorophyll temporal and spatial distribution patterns at basin scales. 3.2 Regional model evaluation for the YS 20 The upwelling is more intense during Spring and weaker in Autumn as indicated in recent observational studies (Ruiz-Castillo et al., 2016). While the model presents upwelling during all the simulated months, this seasonal behavior is well represented in the climatologies obtained from the simulation results shown in Figure 2. The same figure also shows the position of the oceanographic stations occupied during the GOMEX IV cruise, and the delimitation of three areas of particular interest: the inner shelf, the outer shelf, and the upwelling region at Cape Catoche.  ity of the modeled NO 3 is larger than the observed NO 3 at the surface and bottom as shown by the largest standard deviation in Figure 6b.
Below 55 m, both modeled and observed NO 3 are in good agreement in both, the outer shelf and the upwelling area ( Figure   6c). Again, these model results are deemed consistent with observations (since there is no data assimilation) and in the range of other values reported in the literature (Merino, 1997). Since our knowledge, this is the first modeled study focus on the nitrate

Total Nitrogen budget and cross-shelf transports in the YS
The spatial averages of TN suggest a positive trend over the nine simulated years. This is observed for both the inner and the ). This is related to a very efficient biological cycle in terms of TN in the inner shelf, and to the fact that the integration in the outer shelf includes a large volume below the euphotic zone. The temporal series of the integrated 20 TN show a combination of low frequency variability associated with the seasonal cycle and a high-frequency signal (periods lower than 30 days).
To understand the high variability in the YS-TN temporal series and to elucidate the contributions from different sources, the quantification of the cross-shelf fluxes becomes necessary. Their impact on the TN budget and the physical mechanisms modulating such fluxes are investigated next.

25
The cross-shelf fluxes are quantified for the two compartments, the inner and outer shelves (Figure 1b), and for all of the boundaries of each compartment. An schematic view of the main incoming and outgoing pathways of cross-shelf TN fluxes is shown in Figure 9. The yearly averages of the spatially integrated cross-shelf fluxes are shown in Table 1. Denitrification is a form of anaerobic microbial respiration in which nitrate and nitrite are finally reduced to dinitrogen. It represents a major sink for bioavailable nitrogen. Our results suggest that denitrification removes up to the 53% of the TN into in the inner shelf, a significant percentage that agrees with estimates for other shelves in the GoM (Fennel et al., 2006). On the other hand, denitrification in the outer shelf only removes 9% of the TN. Our results also indicate that the denitrification rate 25 tends to increase with time for both inner and outer shelves, similar to TN concentration (Figures 8a and c). This is expected since denitrification is a reduction process, hence an increases in the nitrate concentration means more available DIN to be reduced to dinitrogen.

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 30 and organic matter. We suggest that at least two process are responsible for such modulation: the CTWs and the interaction of the Yucatan Current with the shelf break.  1977). The CTWs have a signature in the sea level, that is well captured in relatively high resolution models such as the one used in the present study (∼ 5 km). These waves have been reported as being responsible for the modulation of upwelling systems such as the Australian coastal upwelling (Shaeffer et al., 2014). Few modelling and observational studies (Kolodziejczyk et al., 2011;Dubranna et al., 2011;Jouanno et al., 2016) describe the characteristics of CTWs in the GoM. In this study, the presence of CTWs is corroborated and its effect modulating the cross-shelf nutrient fluxes at the west margin of the YS is exposed.  (Figure 11b).

15
The multi-taper spectral analysis of Figure 12 (left hand figures) shows that the spectral shape and the range between 10 and 13 days of significant period peaks, are in close agreement between SLA and cross-shelf velocity time series (Figures   12a and b). Although similar peak periods between 7 and 12 days are found for the wind-stress time series (Figure 12c), the resemblance with the cross-shelf velocity is not so evident as for the SLA. Moreover, analysis of the wavelet power spectra of the cross-shelf velocity ( Figure 12, right hand figures), sea level anomaly (Figure 12b) and the wind-stress (Figure 12c), 20 shows that characteristic periods of the cross-shelf velocity in the 50 m isobath are well correlated with characteristic mean periods shown by the sea level anomaly temporal series (r 2 =0.41). By contrast, these periods are not correlated with the mean periods exhibited by the local wind-stress (r 2 =0.16). This shows that the ocean perturbations are remotely generated, with small influence from the local wind variability. The characteristic periods of each wavelet power spectrum are shown in Figure   12d. These periods are computed from the spectral wavelet analysis, by connecting periods of maximum power energy per day. 25 The correlation coefficients found are statistically significant with p-values that greatly exceed the 95% confidence level. The mean and standard deviation of each time-series shown in Figure 12d are displayed in the wavelet power spectra. Notice that a mean period of 10.4 days for the cross-shelf velocity wavelet spectrum matches the mean period found for the SLAs. (rather than in the western Pacific as the classic El Niño) and is surrounded by colder waters to the east and west. As a result wet conditions are produced in the central Pacific. There are evidences suggesting that El Niño Modoki can lead to an increase in the activity of tropical cyclones (Larson et al., 2012). Notice however that the year 2005, widely recognized to be the most intense hurricane season of the GoM, is not marked in the El Niño Modoki index and in the wavelet analysis presented here.
Hence, our analysis is not enough to draw any conclusion and further modeling and sampling effort research is needed.  The interaction of the Yucatan current with the shelf break has been proposed as a mechanism to uplift deep waters into the shelf, feeding the upwelling process (Merino, 1997;Enriquez et al., 2013;Jouanno et al., 2018). Correlation analysis between the Yucatan Current and the TN, PON and DIN fluxes at the eastern margin show high values at interannual time scales ( Figure   13). The trajectory of the Yucatan Current and its closeness to the YS has been also proposed as an upwelling mechanism (Enriquez et al., 2010(Enriquez et al., , 2013Jouanno et al., 2018). When the Yucatan Current is closest to the Yucatan Peninsula, an increase 25 in the sea surface height can be produced, favoring that Caribbean waters enters into the shelf. We have computed the closeness of the Yucatan Current core to the Cape Catoche and Notch areas, which are the two places where nutrients can be upwelled (not shown) (Jouanno et al., 2018). Our results pointed out that there is not a seasonality in the closeness of the current that explains the intensity of the upwelling during Spring or its weakness during Autumn. However, we found an smaller but significant correlation with the sea surface height, supporting the reasoning of a dynamical upwelling due water accumulation 30 in the coast. This analysis also show that the slope of the isotherm 22.5°C, which traces the upwelled water (Cochrane, 1968;Merino, 1997), is slightly more steepness when the current core is close to the shelf.
This study suggests that, at the YS eastern boundary, the nutrient flux towards the coast can be driven by the Ekman bottom layer transport from the interaction of the Yucatan Current with the shelf break. The stress exerted by the intense along-shore velocity of the Yucatan Current at the sharp shelf break produce a bottom generates an Ekman spiral at the bottom. Similar to a surface Ekman spiral in the Northern Hemisphere, a northward along-shelf flow will generate a net depth integrated transport to the left in the bottom boundary layer, i.e., a cross-shelf transport towards the shelf. Similar mechanism has been for the southeastern Australian shelf Shaeffer et al. (2014). Using glider observations, the authors found that the bottom Ekman transport can explain up to the 71% of the bottom cross-shelf transport variability. It is worth noting that Sheinbaum et al.

5
(2016), throughout three years of mooring measurements, show that the mean current near 1000 m oriented along the eastern edge of the YS flows towards the Caribbean sea, i.e., southwards. A southward flow at the eastern part of the YS will generates a coastal downwelling by an eastward Ekman transport. Notwithstanding, they also found that near surface, between 60 and 500 m, the mean current direction is northwest and its magnitude is stronger than the deepest southwards current.
Until now this process has not been fully tested. The present study provides a first modeling evidence that a bottom Ekman 10 mechanism can be the precursor for the upwelling in Cape Catoche. The Bottom Ekman Transport (U bE , m 2 s −1 ) can be and v b are the bottom velocities at the 250 m isobath, f the Coriolis frequency and ρ o = 2015 kg m −3 the reference potential density of the sea water. We found that U bE is positive towards the shelf, and it is largely correlated with the bottom cross-shelf water transport (r 2 = 0.71, ci = 95%) flowing towards the shelf (Figure 14a). The correlation is also 15 large along time (r 2 = 0.78, ci = 95%) as shown in Figure 14b. The comparison between the bottom layer Ekman transport and the time mean vertically integrated TN transport across the eastern 250 m isobath (Figure 14c) show that the Ekman transport is responsible for 65 % of the TN that is entering the shelf. The remaining flux in upper layers seems must be related to recirculation of the Yucatan Current by baroclinic meanders and small-scale structures whose role on the upwelling process needs to be further investigated. 20 Despite being an important mechanism acting to uplift nutrient-rich waters to the shelf helping the Cape Catoche upwelling, bottom Ekman transport is not the only possible mechanism. Other processes, such as the interaction of the Loop Current with the notch could develop intense vertical velocities thus enhancing the upwelling (Jouanno et al., 2018), encroachment of the current jet and mesoscale structures effects Middleton, 2002, 2004). Regarding current-topography interactions, Jouanno et al. (2018), by running a high-resolution model of the GoM, found that removing the notch from the bathymetry, the 25 upwelling could be reduced in a 50% .
Here, we suggest that the upwelling of Cape Catoche is maintained by the combined effect of physical processes of topographic nature, in which the Ekman bottom transport has an important role, but is not the unique. Other processes, such as, current encroachment and eddy induced upwelling can have also important roles that must be analyzed in future research. 30 We present the results of a nine year simulation from a physical-biogeochemical coupled model for the GoM, focusing on the YS. The TN budget, main nutrient transport pathways and their modulation by physical process over the Yucatan shelf are evaluated. Our work provides a first general view of the shelf physical-biogeochemical coupled system, schematized in Figure   15.

Concluding remarks
The results indicate that TN, especially DIN, enters the outer-shelf through the southern and eastern margins. The TN is then driven by a westward shelf current and is exported to the deep GoM and Campeche Sea through the northern and western boundaries, respectively. In the inner-shelf, the biogeochemical nitrogen-based cycle seems to be effective for NO 3 5 remineralization/consumption by the phytoplankton converting 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 Yucatan peninsula. The denitrification represents the main sink of nutrients for the inner shelf, removing more than the 50% of the nitrogen. Despite the inner shelf contributes to the TN to the west boundary of the outer shelf, this contribution is less than 2%, indicating low connectivity from the inner to the outer Yucatan shelf. On contrary, the outer shelf acts as the main nitrogen supply to the inner Initial estimates carried out here suggest that the bottom Ekman layer transport explains the deep TN flux through the eastern YS boundary since there is a positive mean transport (into the shelf) over the nine simulated years along the the eastern shelf break. This indicates that the friction generated between the Yucatan Current and the shelf break can produce a bottom Ekman spiral with a net transport towards the shelf. Bottom Ekman transport can be arrested by stratification and may not be dominant 25 everywhere along the YS east coast as has been documented in other western boundary upwelling regions (e.g., Middleton, 2002, 2004). This study is focused in the Yucatan shelf region whose hydrodynamics and biogeochemical outputs are previously validated.
However, a general evaluation of the physical model ensures that the main dynamics of the whole GoM is correctly represented.
Thereby, the physical model is also evaluated in terms of their hydrography and general dynamics. Temporal averaged maps 5 of eddy kinetic energy (EKE, m 2 s −2 ) obtained from AVISO sea surface height are used to compare with the EKE computed by the model (Figures A1a and b) for 17 years (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012). The model is able to capture the eddy field exhibited by the altimeter product. The main structure of the Loop Current is well represented by the model. In particular mean and variability of eddy kinetic energy field is reasonably captured by the model. A recognized feature is the hook-like pattern of EKE in the western part of the GoM, between 24 and 28°N (Gough et al., 2019), and is more evident in the standard deviation of the 10 EKE model (Figures A1c and d). This pattern is verified in both simulations and observations and is associated with a strong anticyclonic western boundary current that isolates the western continental shelf from the open ocean GoM. On the other hand, the enhancement of the EKE magnitude in ROMS, particularly at the Yucatan channel and Florida strait, is due to the higher resolution of the model (∼5 km) respect to the altimeter product (∼28 km).
Seasonal climatologies of the sea surface temperature (SST) are also compared with satellital product of Aqua-MODIS 15 atmosphere products (http:// modis-atmos.gsfc.nasa.gov) ( Figure A2). The model SST shows a good agreement with both the interannual and seasonal cycles exhibited by the satellite. 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 coast. The model tends to underestimate the coastal SST during winter and spring times, while overestimates it during summer and autumn, respect to the satellital data. Nevertheless, these differences are less than 0.5°C, and its average is in the order of 0.05°C with an standard deviation of 0.4°C. One of the 20 reasons for this discrepancy lies in that satellite-derived provided the temperature relative only to a few microns of the sea surface, which is in fact named as the skin layer temperature (McKiver et al., 2016). Other reason can be also due to the representation of rivers, specially in the southern part of the GoM, where there exist a lack of temporal observations of the main freshwater sources. The river runoff can modify the surface thermal stratification inducing differences between the skin SST and the observed SST below 0.5 m approximately (Donlon et al., 2002). 25 In order to evaluate the mixed layer depth a total of 2629 ARGO floats profiles, availables in the time period between 1995-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 behavior, since the Gulf is an oligrotrophic region in which primary production is controlled by the vertical advection of nutrients 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 dark nutrient-rich waters and the upper ocean layer where the availability of 30 light promotes the growth of phytoplankton and hence zooplankton. Figure A3 shows that the model is able to 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 also shows shallower mixed layer depths during summer, and deepens during winter times respect to the Argo floats. The higher variability of the observed data are likely related with mesoscale structures 13 Biogeosciences Discuss., https://doi.org/10.5194/bg-2019-106 Manuscript under review for journal Biogeosciences Discussion started: 26 April 2019 c Author(s) 2019. CC BY 4.0 License. and submesoscale process which can locally deepening/shallowing the depth of the mixed layer (e.g., Boccaletti et al., 2007;Fox-Kemper et al., 2008;Levy et al., 2012). Moreover, physical processes developed below the 5 km of horizontal resolution are barely permitted by this model configuration. Despite the differences found, the bias between observations and model mixed layer depths are in the order of 1.4 m.