Offshore Freshwater Pathways in the Northern Gulf of Mexico: Impacts of Modeling Choices

The Gulf of Mexico is a very productive and economically important system where riverine runoff acts as a linkage between the continental shelf and the open ocean, providing nutrients in addition to freshwater. This work investigates the three-dimensional transport and pathway structure of this river runoff offshore the continental shelf using ensembles of numerical simulations with different configurations regarding grid resolution (mesoscale resolving and submesoscale permitting) and river setup using suites of 5-months long integrations covering nearly 3 years. The riverine forcing is applied only at the surface over an area around the river mouth, a strategy often adopted in numerical studies, or as a meridional flux with a vertical extension. The simulated flow captures the southward offshore transport of river runoff driven by its interaction with the largest mesoscale circulations in the basin, the Loop Current and Loop Current eddies. This pathway is strong and well-document during summer but also active and relevant in winter, despite a less obvious surface signature. The most intense transport occurs primarily at the peripheries of the Loop Current and the detached eddies, and the freshwater is subducted as deep as 600 m around the mesoscale anticyclonic eddies. Submesoscale motions strengthen slightly the spread of freshwater plumes in summer but their contribution is negligible, if not negative, in winter. Differences in the freshwater distribution and transport volume among runs are small and generally less than 10% among ensembles, with overall slightly higher volume of freshwater transported off-shore and at depth in submesoscale permitting runs that include a velocity flux in their riverine input representation.


INTRODUCTION
Understanding the pathway and transport of freshwater input into an ocean basin is key for estimating land-ocean linkages, the contribution of coastal areas to the oceanic carbon cycling and one component of the anthropogenic pressure on marine ecosystems. In coastal areas, riverine inflow influences oceanic physical and biogeochemical properties. On the physical side, the inflow of low salinity waters modifies stratification, stability and at times circulation (e.g., in the South China Sea, Cai and Gan, 2017;Zeng et al., 2021), while on the biogeochemical side algal blooms, hypoxia, eutrophication and pollutants accumulation result from nutrient and sediment loading by the riverine water, often modified by anthropogenic stressors (e.g., Krug, 1993).
The Gulf of Mexico (GoM, Figure 1A) is a semi-enclosed region highly impacted by several large rivers and is an ideal testbed for studies on river discharge, plume formation, and freshwater transport. The Mississippi River and the neighboring Atchafalaya River (MARS, Figure 1B) are the first and secondlargest freshwater inputs into the GoM, with an average rate of more than 14,000 m 3 /s, and constitute the biggest river system in North America (Morey et al., 2003;Hu et al., 2005;Schiller et al., 2011). They contribute freshwater and abundant nutrients to the semienclosed basin, and this input also enriches microbial diversity (Mason et al., 2016). Rivers smaller than the MARS also drain in the GoM, and their contribution cannot be neglected because of their specific input locations or seasonality. For example, the Mobile River System (confluence of Alabama and Tombigbee rivers) and the Texas Gulf Coast Hydrologic Region (Sabine, Brazos, Trinity, etc.) can contribute to the GoM as much as 3,500 m 3 /s and 1,400 m 3 /s discharge in March and May-June, respectively (Fournier et al., 2019).
Many studies have shown that the advection of riverine water offshore into the northern GoM is modulated by the seasonally varying wind patterns (Walker et al., 1994;Schiller et al., 2011) and by entrainment events due to occasional intrusions of the Loop Current (LC) into De Soto Canyon and the accompanying LC eddies (Morey et al., 2003;Hu et al., 2005;Schiller and Kourafalou, 2014;Androulidakis et al., 2019;Brokaw et al., 2019). These studies are predominantly focused on the interaction of the river plume and the local mesoscale LC eddies during late spring and summer, when the riverine contribution is generally at its maximum, and the impact of submesoscale currents has not been explored in detail.
Global studies have investigated how representing rivering fluxes as bogus volume fluxes imposed at the ocean surface or point sources may impact global salinity budgets (Feng et al., 2021) but less is known for enclosed basins, where models can be run at high resolution to better represent the dynamics responsible for the offshore and vertical spreading of the freshwater anomalies. Among these dynamics, in the past two decades the oceanographic community has paid particular attention to submesoscale circulations. In the GoM, for example, mesoscale structures and rivering inputs have been shown to jointly contribute to the juxstaposition of water masses with different densities and therefore to the development of frontal and baroclinic instabilities. The end result is the formation of submesoscale lateral convergence zones (Luo et al., 2016;Barkan et al., 2017a,b) where fresh water, nutrients, and pollutants accumulate (Huntley et al., 2015;Choi et al., 2017;Bracco et al., 2019;Justic et al., 2021).
Here, we focus on the representation of the freshwater contribution in an ocean model aiming to answer two questions in the context of regional submesoscale permitting and mesoscale resolving high resolution simulations: (1) Do submesoscale circulations contribute to the pathway and transport characteristics of the freshwater plume in conjunction with the LC system (LCS), and if so how much?

and
(2) Is there a significant difference in the model representation of this pathway whenever the freshwater input is imposed as a volume flux at the river mouths or more simply as an additional passive contribution at the ocean surface (i.e., is treated similarly to a precipitation term)? We stress that some ocean models use the second, simpler strategy, also when run for prediction and predictability studies.
To answer them, we use ensembles of ocean model runs at two horizontal grid resolutions and with different river configurations, and we examine the three-dimensional structure of river discharges from ten major rivers in the GoM ( Figure 1A). The integrations cover the period December 2013 -August 2016, allowing us to explore the variability of the MARS-LCS interactions under different conditions ( Figure 1B). This time interval covers a year of relative low discharge in 2014 (medium winter and low summer discharges), a year of typical discharge in 2015 (low winter, high late spring/early summer discharge), and a year characterized by extreme flooding conditions in January-February due to the influence of the El Niño event in 2016 (very high winter and spring discharge).

MATERIALS AND METHODS
To solve the fine-scale, regional circulation, we use the primitiveequation, free-surface, split-explicit Regional Ocean Modeling System (ROMS) in its Coastal and Regional Ocean COmmunity model version (CROCO; McWilliams, 2003, 2011;Debreu et al., 2008Debreu et al., , 2012. The model domain covers the GoM north of 24 • N ( Figure 1A), with open boundaries to the south and east. K-profile parameterization (KPP) is used for vertical mixing of tracers and momentum (Large et al., 1994). The 3rd-order upstream biased advection scheme and the split and rotated 3rd-order upstream biased advection scheme are used for lateral momentum and tracer advection, respectively. Analysis products of the Hybrid Coordinate Ocean Model-Navy Coupled Ocean Data Assimilation (HYCOM-NCODA) experiment 31.0 (before April 4th, 2014) and experiment 32.5 (after April 4th, 2014) at 1/25 • resolution are used as initial and boundary conditions (Cummings, 2005;Cummings and Smedstad, 2013). The boundary conditions are updated every 3 h. All forcing fields (i.e., 10 m wind, surface shortwave radiation, surface downward long wave radiation, precipitation, 2 m air temperature and relative humidity) are from the NAVy Global Environmental Model (NAVGEM), which are also used in HYCOM-NCODA analysis. Atmospheric forcing is employed in our simulations with bulk parameterization, as in HYCOM-NCODA. The topography is derived from 2-min Gridded Global Relief Data (ETOPO2) (Smith and Sandwell, 1997) and is smoothed with a maximum slope factor of 0.25 to reduce horizontal pressure gradient errors caused by sigma-coordinates FIGURE 1 | (A) Bathymetry of the Gulf of Mexico with indicated the river mouth of the rivers considered in the simulations (green dots) and major shelves. Rivers Alabama and Tombigbee share the same river mouth in our model setup. The dots cover the region where the bogus freshwater fluxes are applied in the passive configuration. (B) Freshwater discharge rate of each river from December 2013 to August 2016. Vertical dash lines show the 1st day of January and July during this period. Green dotted line indicates the mean discharge rate, calculated as monthly means over the 11-year period Jan 2010 -Nov 2021. (Sikirić et al., 2009). The configuration choices insure that the Loop Current dynamics do not diverge significantly from the analysis in all runs and that overall the model follows the dynamical evolution of the HYCOM-NCODA system. River discharges of ten major rivers (Figure 1), collected daily by the U.S. Geological Service (USGS) and the U.S. Army Corps of Engineers (USACE), are imposed in the simulations in an active or passive manner. In the active configuration, the discharges are imposed as a southward volume flux at the river mouths. The river salinity is set to 4 psu, and the temperature is set at the local air temperature value. The freshwater flux decays exponentially away from the surface toward the ocean bottom. In the passive river configuration, which is commonly preferred because easier to implement, the river discharge is simulated as a bogus flux that is added to precipitation around the river mouth, with the amplitude of the flux decreasing according to a Gaussian function with an e-folding length of 50 km. (see Figure 1) to ensure that the freshwater reaches realistically far offshore (Barkan et al., 2017a).
The December 2013 -August 2016 period is divided into 8 intervals, each 5 months long. For each interval, the model is run freely (without any data assimilation) starting on December, April and August 1st, 2nd, and 3rd of each year from the HYCOM analysis. The three-member ensembles allow for quantifying the internal variability, while the re-starting every 5 months assures that the runs do not diverge much from the observed conditions in their LCS representation. At the same time, the set-up choices limit us to investigating subseasonal to (at most) seasonal changes associated with the riverine fluxes representation. Two ensembles (active and passive) are run at 1 km (submesoscale-permitting, SP) and two at 3.5 km (mesoscale-resolving, MR) horizontal resolution, with 70 sigma layers in all cases. In the SP case a twomember ensemble is run also without riverine discharge except for that contained in the initial conditions; in those no-river runs the surface salinity continues to be modified by evaporation and precipitation.

RESULTS
The predictability of surface fields -sea surface temperature (SST), salinity (SSS), relative vorticity and height (SSH) -is described in terms of variability, climatology and ensemble spread in Sun et al. (2021), thereafter SBL21, and is not repeated here. In this work, we focus on salinity pattern at and below the surface, paying closer attention to the 3D distribution of riverine waters as they move offshore into the open ocean. A configuration close to that adopted here has been validated in various fields in previous studies (see e.g., Cardona et al., 2016;Luo et al., 2016, andLiu et al., 2021). Overall the model reproduces well the stratification and mesoscale variability, but LC structures (the LC and the detached eddies) are often stronger than observed in terms of their eddy kinetic energy. This is a common problem in ocean-only model. Renault et al. (2016), for example, assessed the role of atmospheric coupling using oceanatmosphere coupled and uncoupled simulations, and showed that FIGURE 2 | Monthly averaged surface salinity in the first member (initialized on the 1st of the month) of the SP active case shown in the third (February, June, and October) and fifth (April, August, and December) months of the (free) simulation. Black contour indicates 0.17 m sea surface height anomaly (SSHa) and broadly follow the LCS and the Rings, following Leben and Born (1993).
the current-wind feedback in coupled simulation attenuates the mean transfer of energy from the atmosphere to the ocean, and consequently decreases the level of eddy kinetic energy in the ocean by nearly 30%.

Modeled Salinity and Circulation
Features in the Gulf of Mexico Figure 2 shows the monthly mean sea surface salinity (SSS) in the 3rd and 5th months of each model period for the first member of the SP active ensemble. The SSS field varies greatly across years because of the freshwater input, the variations in wind pattern as well as the position of the LC and the Rings (Morey et al., 2003;Schiller et al., 2011;Schiller and Kourafalou, 2014;da Silva and Castelao, 2018). Winds in the GoM are typically southereasterly in summer and northerly and northeasterly during fall and winter (Morey et al., 2003) and drive the seasonal circulation over the shallow inner shelf. In the multi-year average, in February -April low-salinity waters are confined to the Louisiana-Texas (LATEX) continental shelf, and extend offshore in summer (June to August), with a major branch of freshwater spreading to the east of the Mississippi Delta, along the coast of the Mississippi-Alabama-Florida (MAFLA) and toward the West Florida shelf. The shape and extension of this branch are linked to the position of the LC. In fall the fresh SSS retracts toward the shelf, with remnants around the LCS, and in winter (November to January) only the shelf areas are occupied by low salinity waters near the surface.
In the period considered, low-salinity waters (identified by salinity less than 34 psu) were mainly limited to the coastal areas in August 2014 because the discharge in April-July was anomalously low, but extended well offshore to as far as 26 • N in 2015, when the late spring/summer discharge was highest among the 3 years ( Figure 1B) and the LC was in an extended position for most of the time.
During winter and spring, the El Niño Southern Oscillation (ENSO) is an important driver of interannual variability in salinity and coastal circulation (Gomez et al., 2019). ENSO indeed modulates the global precipitation patterns over the continental United States that feed the Mississippi River basin and the wind anomalies that control the alongshore currents over the LATEX shelf via the thermal wind relation (Childers et al., 1990;Schimidt et al., 2001). This modulation results in flooding conditions and wider than normal offshore spreading of low-salinity waters in winter and early spring during El Niño events, as verified in 2016 from February to April.
As mentioned, the GoM circulation in the upper 600-800 m of the water column is dominated by the anticyclonic Loop Current (LC), which is relatively warm and salty, and its detached eddies or Rings (Hamilton, 2009). The LC and the Rings can extend to 200 -300 km in diameter (Gordon, 1967). Whenever a Ring is formed, and through its westward drift, smaller mesoscale eddies, spinning both cyclonically and anticyclonically, are also generated. Around and within the large mesoscale structures, submesoscale circulations have been observed year-round (Poje et al., 2014;D'Asaro et al., 2018) at horizontal scales from hundreds of meters to few kilometers. They emerge from the conversion of available potential energy from the mixed layer to kinetic energy and are fueled by the density gradients associated with the energetic mesoscale structures (McWilliams, 2008) and the river discharge (Barkan et al., 2017a). Figure 3 presents the surface salinity gradient field, the surface relative vorticity ζ = ∂v/∂×−∂u/∂y normalized by the Coriolis parameter f, and the near-surface stratification, or the squared Brunt-Väisälä frequency N 2 = − g ρ dρ dz , where ρ is the potential density, in 2015 for one of the SP active runs. Similar patterns are detectd in other cases but with relatively weaker magnitude or less extension for all fields. For relative vorticity, in the offshore portion of the basin -offshore the shallow continental slope -submesoscale circulations in the form of small eddies are most energetic in winter when the mixed layer is deepest and submesoscale fronts prevail in summer, due to the river-induced submesoscale enhancement, as shown in previous studies (Luo et al., 2016;Barkan et al., 2017b). Over the shelf, submesoscale eddies are more commonly observed from late spring to early fall (Horner-Devine et al., 2015;Kobashi and Hetland, 2020). Focusing on stratification, very weakly stratified waters can be seen inside the LC and Rings extending to the shelves in Frontiers in Marine Science | www.frontiersin.org winter, since the nearshore temperature is low throughout the shallow water column and winds strongest at the coast. On the other hand, waters are strongly stratified in summer whenever impacted by the abundant riverine output (Morey et al., 2003).
Patterns of salinity gradients, stratification and vorticity at the mesoscales remain considerably similar across ensembles, as seen in Supplementary Figure 1 for the corresponding quantities in the MR ensembles, and in SBL21 for further comparisons of predictability skills. Figure 4 shows the monthly averaged spatial patterns of the freshwater spreading area and normalized volume, or thickness in 2015, together with the sea surface height anomaly (SSHa) contours across the different ensembles. Again, the freshwater is identified by 34 psu isohaline. Sea SSHa contours of 0.17, 0.4, and 0.6 m outline major anticyclonic features, and the 0.17 m threshold is assumed to identify the boundary of the LC and LC eddies in the GoM following previous literature (Leben and Born, 1993;Donohue et al., 2016;Hamilton et al., 2016). In winter, westward/southward wind-driven coastal currents carry low salinity water over the LATEX shelf. During the summer, two branches of low salinity water can be present and transported eastward in the northern GoM. The first branch represents the portion of pre-existing waters on the LATEX shelf transported along-shore, which is also present in the no-river case. The second branch originates mainly from the MARS and moves eastward onto the MAFLA shelf, the West Florida shelf or is advected into the Gulf interior by the intrusion of the LC and the Rings. Finally, in August parcels of low salinity waters can be seen trapped in the cyclonic recirculation between the large anticyclonic patters (LC and Rings) and they are more relevant in the SP-active case than in the other ensembles. The trapping of low salinity waters between two mesoscale anticyclonic structures has been pointed by Androulidakis et al. (2019) as a secondary mechanism of freshwater removal by the LC eddies. The SSHa patterns are overall similar among different ensembles. Quantification of the differences is presented in SBL21, where it is found that the LC has a statistically significant lower probability to extend northward of 28 • N in the SP integrations with an active river configuration, compared to the other cases. In terms of salinity, the modeled freshwater spatial distribution is modestly affected by the setup of river forcing. This is the case although Morey et al. (2003) have demonstrated that the variability in river discharge does not control, but only modulate, the seasonal salinity patterns in the GoM, with the LCS being a player at least as important. In our simulations, both the spreading area and volume are small in the no-river case, as to be expected, and maximum values are achieved in the active river simulations. The role of small differences in LCS are most evident in comparing the SP passive and active cases and their salinity distribution in August 2015. The extension of the LC, which is slightly different among the two cases, modifies the pathway of the Mississippi River water pushing it further east in the passive case.

Interannual and Seasonal Variability in Freshwater Plume
To quantify the different behaviors in the transport pathways and transport efficiency, we calculated the time-series of freshwater spreading area and volume, shown in Figure 5. We excluded from the calculation coastal areas where the water column depth is less than 100 m given that differences could be driven by river  configuration choices in the immediate surrounding of the river mouths. This arbitrary threshold excludes all areas impacted by the modified precipitation fluxes where waters may stagnate in conditions of low winds. The Mississippi River mouth represents the exception, but its freshwater is rapidly transported offshore or to the west in both passive and active cases. We varied the threshold and tested values between 50 and 200 m and results do not change qualitatively.
Interannual changes are generally larger than the differences among different ensembles. As a reminder, for each 5-months long simulation period, the integrations share the atmospheric forcing and start from similar initial conditions (1-day apart). The freshwater area coverage tends to be low in winter, except for 2016, reaches ∼2-3 × 10 11 m 2 during summer, and drops to less than 10 11 m 2 in November of 2014 and 2015. The freshwater volume follows an analogous evolution with comparable seasonality. The largest area and volume of lowsalinity waters occur in the summer of 2015. Moreover, relatively lower freshwater coverage but comparable to 2015 freshwater volume is found in August 2016, implying that low salinity waters reached deeper than in other summers, because the larger fluxes occurred earlier in the year (Figure 1), when the MLD was deeper. In the 2016 summer the introduction of freshwater in the basin due to extreme precipitation events, which is illustrated by the increasing freshwater area and volume in the no-river case, is also a noticeable occurrence.
Differences among ensembles are small. The MR and SP runs have similar LCS representations and identical river forcing,  and the freshwater spreading evolves similarly over time. In the active river configuration the velocity flux associated with the freshwater input contributes to a slightly greater spread compared to the passive set-up. Differences are larger whenever the freshwater discharge is high, as verified during winter/spring 2016 and more generally in the summer season.
The freshwater spreading is a near-surface process strongly affected by the upper-layer circulation dynamics. During late spring and summer the submesoscale circulations tend to further aid spreading in the SP simulations. In these seasons the prevailing submesoscale features are in the form of density fronts with horizontal advective tendency induced by confluence and convergence (Wang et al., 2021) that contribute to the transport of the riverine water. Indeed, in summer the discrepancy between SP and MR simulations is consistent among years. The lateral convergence zones in the submesoscale frontal structures that abound around the LCS (Zhong et al., 2012;D'Asaro et al., 2018) accumulate and distribute freshwater, contributing to the wider extension of riverine inputs. The summer submesoscale contribution is highlighted comparing the instantaneous distributions of surface salinity and frontogenetic tendency, F, defined as Dt is the total derivative and represents the rate of change ( Da Dt = ∂a ∂t + u ∂a ∂x + v ∂a ∂y ) and Q = (Q1, Q2) = −( ∂u ∂x ∂ρ ∂x + ∂v ∂x ∂ρ ∂y , ∂u ∂y ∂ρ ∂x + ∂v ∂y ∂ρ ∂y ) (Hoskins, 1982; Frontiers in Marine Science | www.frontiersin.org  Figure 6). F, which is two orders of magnitude greater in the SP-active case compared to the MR-passive one, follows closely the freshwater pathways in all cases.
From late fall to early spring, on the other hand, the mixed-layer is deeper and offshore the continental shelf fronts and filaments can be more easily fragmented by submesoscale instabilities that extract energy from the mixed-layer and generate eddies that restratify the upper water column (Luo et al., 2016). The time evolution of F over waters deeper than 100 m is shown in Figure 7, where is also evident that within each given ensemble the internal variability is small even at the submesoscales. Given the conbined annual cycle of F and mixedlayer depth, the prevailing submesoscale circulations converge material at their interior (D'Asaro et al., 2018), stretching the freshwater patches in filaments less effectively than in summer. Furthermore, in these seasons low-salinity waters flow mostly along the LATEX shelve, due to the prevailing northeasterly and northerly winds and more limited river discharge. The submesoscale eddies and fragmented fronts that occupy waters with a depth greater than ∼ 100 m, have negligible impact on the freshwater spread in "normal" winters (2014 and 2015). On the other hand, in 2016, when the river runoff was very abundant and was pushed offshore by the anomalous winds, the submesoscale circulations cause a slight reduction in volume spread compared to the MR case by limiting the mixed-layer extention (Barkan et al., 2017b).

Horizontal and Vertical Transport of Riverine Waters
We next analyze the freshwater transport, T, highlighting the riverine water pathways. We do so by adopting the methodology introduced by Schiller et al. (2011) and applied to the GoM in subsequent works (Schiller and Kourafalou, 2014;Androulidakis et al., 2019;Brokaw et al., 2019). The freshwater transport is defined as T = V × fw × dx × dz where V is the horizontal velocity with its zonal or meridional components, dx is the corresponding space distance, dz is the vertical thickness, and fw is the freshwater anomaly calculated as fw = (S b − S) / S b . Here S b is the background salinity field without the added river discharge (i.e., the mean salinity in no-river cases in this work) and S is the salinity field in each of the active or passive runs without any threshold applied. Since V and fw can be both positive and negative, positive zonal transport can be interpreted as either freshwater moving to the east (both positive) or saltier water moving westward (both negative). In the following we modify fw by setting all the negative values to zero in order to select the transport of freshwater. The difference between setting the negative fw to zero and considering freshwater anomalies of both signs is shown in Supplementary Figure 2, while Figure 2 and Supplementary Figure 3 display monthly salinity fields in one of the SP no-river and SP-active members, respectively. We stress that the mesoscale patterns are similar, but necessarily not identical among the runs (see also SBL21), and that the SP no-river ensemble includes contributions from precipitation anomalies in its freshwater budget.
The spatial patterns of freshwater transport integrated over the upper 100 m of the water column in the last month of the 5-month long runs in 2015 are shown in Figure 8 for the first member of the SP active-river ensemble. Positive values correspond to eastward and northward transport. Anticyclonic (clockwise) freshwater transport greater than 100 m 3 /s is detected at the periphery of large anticyclones in summer, especially near the north, east and south boundaries of the LC eddy (Figure 8, transport vector). A similar pattern is also observed in winter but with reduced transport associated with it. A small cyclonic eddy (Figure 3) with opposite transport can indeed be seen at about 94 • W, 25 • N during February.
Previous studies on the offshore pathways of freshwater transport in the GoM have focused on the interactions between the LCS and MARS in summer (Schiller and Kourafalou, 2014;Androulidakis et al., 2019;Brokaw et al., 2019), following the peak of riverine discharge. Brokaw et al. (2019) identified three main pathways in summer by analyzing the surface salinity transport from 2010 to 2018, two of which are captured here and are realized by the relative positions of the LC and LC eddies with respect to the river plume. According to the first pathway, the freshwater is advected eastward from the MAFLA shelf, then southward and finally exits the GoM from the Florida Straits. The second pathway is characterized by narrow freshwater bands surrounding the LC eddies or the LC extension. This recirculation pattern advects the low-salinity waters offshore to the south, then westward and back north toward its source. The third pathway describes conditions when the Rings and the LC do not reach north enough to interact with the river plumes and no detectable interactions between the LCS and the MARS are found. This condition was relevant to portions of 2010 and 2012, when the freshwater remained distributed in the northern part of the basin with decreasing concentrations moving southward and a nearly homogenous gradient, but in the period considered here this third condition is never observed or simulated. On the other hand, a west-to-east offshore transport characterizes the summer of 2015 when two large Rings occupy the central GoM. The smaller cyclonic circulation found at about 93 • W and 26 • N in August acts as a bridge connecting these two anticyclones and contributes to the lateral transport of nearshore freshwaters from west to east. A similar redistribution pattern is observed also in the other two summers (Supplementary Figures 4, 5).
Our analysis indicates that the influence of the LC system on the offshore transport of freshwater extends to the winter seasons, despite its surface signature being significantly reduced. As mentioned, in winter the surface river discharge is confined to the western coastal regions due to the prevailing northerly and northeasterly winds. The low-salinity waters surrounding the LC and central LC eddies mainly originate from the confluence regions of downcoast inner shelf currents and upcoast outer shelf currents between 26 and 28 • N, west of 96 • W (Walker et al., 1996;Zavala-Hidalgo, 2003). The offshore transport is sustained by the jet associated with the decaying eddy or a pair of eddies usually found near the coast in this season (Zavala-Hidalgo, 2003), and seen, for example, between 93 • and 96 • W in Figure 3. This longdistance transport mechanism of nearshore, river-influenced waters is evident in April 2015 (Figure 8) and occurs also in 2014 and 2016 (Supplementary Figures 4, 5). Another important contribution to winter transport comes from summer coastal waters transported offshore that continue to be advected by the mesoscale field.
To expand our analysis to greater depths and fully encompass the GoM deep chlorophyll maximum located at around 100 m (e.g., Pasqueron de Fommervault et al., 2017) and the whole euphotic layer, we extended the computation of the integrated freshwater transport between 100 and 200 m (Figure 9). This calculation is relevant for the primary productivity of the basin being riverine water the main contributor to its nitrogen and phosphorus balance. Below 100 m the freshwater transport follows the boundary of the Rings and LC, with a pattern similar to that of the upper layer but overall decreased magnitude.
Significantly enhanced values occur also inside/around cyclonic structures that are part of the LCS, which in August 2015 can be seen at about 94 • W and 25 • N in February and 93 • W and 26 • N. Transport patterns and overall intensity are slightly greater in active runs than in passive configurations. This is quantified in Figure 10, where the probability density function of the transport vector is plotted for the month of August in all 3 years and two SP ensembles. Differences are small but robust, and corroborate the volume differences in low salinity waters quantified in Figure 5. Indeed, if the tails of the distributions are similar once the ensemble variability is accounted for, the core indicates greater transport in the active ensemble in the 50-300 m 3 /s in all members and all years. In addition, the largest transport rate identified by the flattest slope is detected in 2015 for both active and passive cases, while in 2014, the overall transport rate stays at the lowest level among all 3 years.
Finally, the LCS controls not only the horizontal but also the vertical distribution of freshwater. In the passive configuration, the strength of the subduction around the LC and the Rings is reduced compared to the active runs. This is exemplified in Figure 11, where transects taken at latitudes ranging from 26 to 28 • N are shown for June 2015. Here the choice of the month is due to the very strong similarity in the LCS in the two members considered. The contour of 1026 kg/m 3 density is used to outline the position of different circulation features. The freshwater can penetrate to depths close to 600 m during summer, approaching the vertical scale of the LCS structures. The most intense vertical penetration is found at the periphery of the Rings. This feature, combined with the fact that meridional currents weaken with depth, contributes to a strong downwelling/subduction of freshwater into the deep ocean. In the northern portion of the basin, where the LC and the Rings do not intrude, the freshwater transport below the mixed layer is modest or null.

CONCLUSION
In this work we investigated the three-dimensional pathway and offshore transport of low-salinity riverine outflow in the GoM over a nearly 3-year period, using regional ocean simulations at mesoscale resolving and submesoscale permitting resolution. We extended previous studies that focused on the offshore removal process of near-surface freshwater from the MARS by exploring the intraseasonal variability of the offshore water transport. We considered simulations with different river configurations, ranging from no river inflow (no-river), to a passive configuration in which the riverine fluxes are converted to precipitation and are applied to an area extending offshore about 50 km from the river mouth, and finally to an active configuration with a point-source southward volume flux added at the mouth of each river. The passive configuration is usually preferred and has been used in several recent works showing realistic offshore extension of lowsalinity waters in the GoM (Barkan et al., 2017a;Liu et al., 2021). However, the active configuration proves that converting the river discharge to a vertically decaying volume flux can increase the amount of riverine water that reaches offshore in the model runs, and its vertical extent.
It is well documented that river discharge that originated from the MARS interacts with the LC system in the summer season (Morey et al., 2003;Schiller et al., 2011;Schiller and Kourafalou, 2014;Androulidakis et al., 2019;Brokaw et al., 2019). During summer, the intrusions of the anticyclonic LC eddies entrain at their periphery freshwater from the river-influenced regions and move it toward the Gulf interior. These low-salinity waters are then transported either eastward by the LC and can exit through the Florida Straits, or westward to the central Gulf by the LC eddies. We showed that in winter the lateral transport of nearshore waters from confluence regions over the western shelf is the dominant pathway, but the southward transport of MARS waters by the LC eddies remains significant while less notable since the mixed-layer is deeper, the surface signature consequently less evident, and the amount of involved freshwater generally smaller (Figure 7).
Our analysis also provides evidence that the freshwater transport surrounding the LC and LC eddies penetrates as deep as ∼ 600 m. The low-salinity subduction processes are most prominent around the Rings and in the cyclonic circulations commonly generated when the LC eddies detach from the main current (Figures 8, 10).
On interannual timescales, the rainfall over the continental United States and the LC dynamics modulate the offshore reach of the riverine input in spring to fall, while the ENSOinduced co-variability of salinity and coastal currents in the northern GoM (Gomez et al., 2019) is the major player in winter (Munoz and Dee, 2017). Particularly, El Niño events are linked to large positive precipitation anomalies in late fall and winter, and declining values from spring onward (Lee et al., 2014(Lee et al., , 2018. This rainfall anomaly and the associated storminess tendency that co-occur with strong El Niños, contributed to the extensive three-dimensional spread of low-salinity waters in the winter of 2016. The waters of riverine origin were more abundant and also transported to greater depths in spring 2016 compared to the spring in the other 2 years. In addition to the dominant mesoscale LC circulations, we discussed the role of submesoscale motions in modulating the transport and distribution of freshwater in the GoM. Generally, the intensified submesoscale fronts accumulate, stretch and advect the freshwater plume (Zhong et al., 2012;D'Asaro et al., 2018). Given their abundance in the summer season, higher values of freshwater spreading area and volume are found in the SP cases in summer, following the maximum river runoff, when both mesoscale and submesoscale currents have the strongest influence on the distribution of low-salinity waters. The difference, however, is contained between 3 and 14% in area, and between 0 and 10% in volume for the offshore portion of the domain. In 2016, when the riverine inflow was large also in winter, the restratification by submesoscale eddies had a small but negative impact on the freshwater spreading, so that the overall freshwater volume found offshore is greater in the mesoscale resolving run.
In the vertical direction, submesoscale dynamics can significantly increase the vertical velocity and mixing rate within the mixed layers (Thomas and Ferrari, 2008;Zhong and Bracco, 2013;Callies et al., 2015;Liu et al., 2018Liu et al., , 2021, and accelerate the subduction process of low-salinity waters as a result. Their offshore contribution, however, is limited in summer by the mixed-layer depth . As a result, the subduction in this season is mostly driven by the mesoscale circulations with comparable efficiency in the MR and SP ensembles. Overall, differences are robust but small, with the submesoscale permitting resolution and the active river configuration contributing to the offshore spreading by making it slightly more efficient and deep reaching in most seasons. On the other hand, the predictability of offshore spreading of riverine water in all years considered is intrinsically linked to the internal (modeled) variability of the LCS and does not change significantly among ensembles. This is good news for biogeochemistry studies. While submesoscale circulations may impact the distribution and availability of nutrients, their contribution to the offshore redistribution of riverine waters remains small even in presence of very large outflows as in the case of the MARS.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.dropbox. com/sh/g5sywqcyqog0i8v/AAD2LvGMlCy2rHk-t1zQzHA8a? dl=0.

AUTHOR CONTRIBUTIONS
AB conceived the work. DS performed the simulations. GL conducted the analyses and wrote the first draft. All authors contributed to the interpretation of the results and to the writing of the article, to the article and approved the submitted version.

FUNDING
This work was supported by the Gulf Research Program (GRP) of the National Academies of Sciences, Engineering, and Medicine through its Understanding Gulf Ocean Systems 2 (UGOS-2) Grant "The Loop Current and the Mississippi-Atchafalaya River System: Interactions, Variability, and Modeling Requirements." GL and AB were partially supported by the National Science Foundation (OCE -1658174).