Simulating Water Residence Time in the Coastal Ocean: A Global Perspective

Exchanges between coastal and oceanic waters shape both coastal ecosystem processes and signatures that they impart on global biogeochemical cycles. The timescales of these exchanges, however, are poorly represented in current‐generation, coarse‐grid climate models. Here we provide a novel global perspective on coastal residence time (CRT) and its spatio‐temporal variability using a new age tracer implemented in global ocean models. Simulated CRTs range widely from several days in narrow boundary currents to multiple years on broader shelves and in semi‐enclosed seas, in agreement with available observations. Overall, CRT is better characterized in high‐resolution models (1/8° and 1/4°) than in the coarser (1° and 1/2°) versions. This is in large part because coastal and open ocean grid cells are more directly connected in coarse models, prone to erroneous coastal flushing and an underestimated CRT. Additionally, we find that geometric enclosure of a coastal system places an important constraint on CRT.


Introduction
Coastal seas differ from the open ocean in numerous ways, tending to be disproportionally shallow and strongly influenced by benthic and terrestrial processes. Their contributions to total global primary production (20-30%; Wollast, 1991;Longhurst, 1995), organic matter burial (80%), sedimentary remineralization (90%; Gattuso et al., 1998), and calcium carbonate deposition (50%; Pernetta and Milliman, 1995) far outweigh their relatively modest contribution (7-8%) to global oceanic area (de Haas et al., 2002). These shallow and productive waters host diverse and strongly interactive pelagic and benthic ecosystems, serve as important habitats for juvenile fish recruitment (Mieszkowska et al., 2014;Miloslavich et al., 2011), and support the vast majority of global fish catch (Pauly & Zeller, 2016;Stock et al., 2017). Complex hydrodynamics in the coastal ocean shape not only biogeochemical processes within it but transport pathways of nutrients and carbon from land to the deep ocean (Bauer et al., 2013;Bauer & Druffel, 1998;Cai, 2011).
Anthropogenic activities across the globe have resulted in substantial changes that threaten the functionality of coastal habitats (e.g., Diaz & Rosenberg, 2008;Heisler et al., 2008). Over the past decades, climate models have been used as powerful tools to explore anthropogenic impacts on marine ecosystems at regional to global scales (Bopp et al., 2013;Stock et al., 2011). However, with coarse horizontal grids (1-3°) these models are invariably unable to resolve the complex coastal bathymetry and fine-scale physical processes relevant to nutrient cycling, primary productivity, and other key biogeochemical processes in the coastal ocean (Holt et al., 2017). Thus, detailed studies on coast-ocean exchanges, the role of coastal processes in global biogeochemical cycles, and their responses to climate change have been mostly limited to regional studies (e.g., Glibert et al., 2014;Hermann et al., 2016;Holt et al., 2015;Laurent et al., 2017).
A fundamental property of a coastal system is the amount of time its waters stay in that environment. This retention timescale constrains coast-ocean exchanges and allows for biogeochemical processes to impart their unique signatures on coastal environment. It has been shown that coastal retention may enhance productivity (Edman et al., 2018;Lohrenz et al., 1999), favor recruitment of larval fish (Mullaney & Suthers, 2013), constrain connectivity between coral reefs (Figueiredo et al., 2014), and facilitate development of harmful algal blooms (Pitcher et al., 2010) and coastal hypoxia Rabalais et al., 2014). Coastal retention and coast-ocean exchanges are also key to the role that the coastal ocean plays in the global carbon cycle. For instance, Bourgeois et al. (2016) suggests that coastal sinks of anthropogenic CO 2 are weaker (~4%) than expected from the area because the slow exchanges are inadequate to transport the excessive CO 2 from coastal waters to the deep open ocean.

10.1029/2019GL085097
• Table S1 • Movie S1 Correspondence to: X. Liu, xiaoliu@princeton.edu Water residence time has been previously investigated using both observational (Colbert & Hammond, 2007;Moore, 2007) and numerical (Delhez et al., 2004;Shen & Haas, 2004;Zhang et al., 2010) approaches. While earlier studies on residence time are mostly limited to regional scales, several global syntheses have been made recently. For example, Bourgeois et al. (2016) used a 1/2°global ocean model to compute residence time for each considered coastal domain by dividing its volume by the integrated outflow at the coast-ocean boundary. Sharples et al. (2017) uses hydrodynamic principles to derive residence times of freshwater plumes on the shelves as a function of latitude and the shelf width.
Here we proposed a novel, tracer-based approach to quantify coastal residence time (CRT) globally. Using a suite of newly developed global ocean-ice models (at 1°, 1/2°, 1/4°, and 1/8°), this work provided a global perspective on model simulated CRT and its spatio-temporal variability. We also assessed the sensitivity of simulated CRT to model resolution and explored the importance of geometric and geostrophic constraints on CRT. We focused our analyses on the coastal domains of 60 large marine ecosystems (hereafter referred as the coastally restricted LMEs or cLMEs; http://lme.edc.uri.edu) that are hotspots of nutrient overenrichment, marine productivity, and fisheries catch (Conti & Scardi, 2010;Friedland et al., 2012;Stock et al., 2017).

The Global Ocean-Ice Models
The coupled, global ocean-ice models were configured from the latest generation of the Modular Ocean Model (MOM6) and the accompanying Sea Ice Simulator (SIS2) both developed at the NOAA GFDL (Adcroft et al., 2019). The model suite includes four versions of MOM6-SIS2 with a nominal horizontal resolution of 1°(~100 km), 1/2°(~50 km), 1/4°(~25 km), and 1/8°(~12.5 km), respectively (see supporting information Text S1 for further details). All four models were initialized at rest using long-term climatologies of observed temperature and salinity (World Ocean Atlas 2013 v2). Experiments were run for a total of 40 years with two repeating cycles between 1988 and 2007 following the CORE forcing protocol (Large & Yeager, 2009). We considered the first 30 years as the spin-up period to allow coastal circulation and tracers to reach equilibrium and only analyzed model output from the last 10 years (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007).

An Idealized Tracer to Estimate Coastal Residence Time
Here "coastal residence time" (CRT) is defined as the elapsed time since a parcel of source water enters the coastal domain (defined by the 200-m isobath) before it exits to the open ocean while allowing excursions across the coast-ocean boundary. A novel aspect of this definition is that the global coastal ocean is viewed as a connected water body such that CRT records the total time a water parcel stays in any part of the global coastal ocean rather than a specified domain (i.e., a water parcel would accumulate CRT while traveling alongshore from one coastal system to another). To estimate CRT, we implemented an idealized coastalage tracer in the models. This method represents a novel, coastal application of an established tracer approach commonly used to estimate ventilation timescales of water masses (England, 1995). Specifically, at release, the age tracer concentration (C age ) starts increasing with time in all grid cells within the coastal domain (C age, t + Δt = C age, t + Δt) and decreasing to a minimum of zero in all grid cells within the open ocean domain (C age, t + Δt = min(0, C age, t − Δt)). The purpose of allowing a gradual decrease of Δt in the open ocean, instead of resetting C age to zero, was to prevent fine-scale spatio-temporal variability at the coastocean boundary from degrading the representativenss of C age for coastal retention. As time elapses, C age at a given coastal location will approach equilibrium when the gain due to "aging" is balanced by advective and diffusive losses to the open ocean. At equilibrium, C age is on average lower at locations near the open ocean boundaries, where coastal waters are exchanged rapidly with the "younger" (lower C age ) oceanic waters, and generally increases toward the coastal interior, where waters interact with the "older" (higher C age ) along-shore currents. The equilibrium concentration of C age is thus a measure of the amount of time an average water parcel in that given location has remained in the coastal domain or, in other words, its CRT.
To assess the spatio-temporal variability of CRT, we focused on cLMEs that are geographically and ecologically linked domains along the global coast. Of the 66 commonly defined systems, we excluded six cLMEs (#10, 23, 62, 61, 63, 64) that are inland seas, island chains, or heavily ice-covered systems. We calculated the characteristic CRT for each remaining 60 cLMEs as the depth-and area-weighted logarithmic mean of C age across all coastal grid cells within that system. The use of logarithmic means ensures that the characteristic CRT for a cLME similarly weights low (high) residence time grid cells in the exterior (interior) of the domain.

Evaluation Datasets 2.3.1. Observed CRTs
Previous studies characterizing water retention time have been done for a limited number of coastal and nearshore systems, either from direct observations on currents and water properties, or based on mass balances of measured geochemical tracers (Table 1 and references listed within). These studies allowed for a first-order assessment of our simulated CRTs (i.e., whether model estimates and observations are of the same order of magnitude). However, a global-scale evaluation was impossible from these limited studies alone. An alternative, globally applicable metric based on the maintenance of contrasting salinity fields between fresher coastal waters and saltier oceanic waters is thus applied here (Session 2.3.2).

Observed Salinity Fields and Freshwater Plumes
Salinity is not a direct measure of flushing rates and residence time but has been successfully used as a proxy for understanding lateral transport in estuaries and coastal seas (e.g., Halverson & Pawlowicz, 2008;Mudge et al., 2008). For example, insofar as freshwater supply to the coastal ocean is adequately represented in a model, a "salty" bias in coastal waters is indicative of the model overestimating coast-ocean exchanges and underestimating CRT, and vice versa. Here we compiled all surface (0-15 m averaged) salinity measurements from stations, CTD, gliders, and floats between 1988 and 2017 from the World Ocean Database (n = 2,179,491; Figure S1). These data were median binned on a 1°×1°global grid seasonally to minimize biases associated with duplication and clustering. We then assessed the observed and modeled surface salinity contrasts between the offshore and coastal domains of each LME separated by the 200-m isobath (calculated as SS offshore − SS coast and hereafter referred to as "coastal salinity deficit"). Additionally, freshwater plume structures were examined for the Northwest Atlantic shelf and the Gulf of Mexico where available observations are relatively rich.

Evaluation of Simulated Coastal Residence Times (CRTs)
All four models predict a wide spectrum of CRTs for the 60 cLMEs globally (Figure 1). In the highestresolution version (1/8°), for example, these CRTs range from 7 days within the Somalia Coastal Current to over 5 years in the East Siberian Sea, with a median of 93 days and an area-weighted, logarithmic mean of 116 days (Figure 1; Table S1; Figure S2). Models are overall in good agreement with observations (Table 1). For example, a simulated CRT of 8-17 days for the California Current suggests that coastal waters in this system are exchanged vigorously with the oceanic waters, consistent with residence time inferred from radium isotope measurements (≤18 days). The Laptev Sea, in contrast, has a much longer CRT, as revealed by estimates from both observations (3-6 years) and models (4-5 years). In general, the 1/8°and 1/4°models perform better in simulating the observed CRTs, in the sense that model estimates for all nine regions differ by less than 50% from observations. The sign and magnitude of model-observation differences vary both by cLMEs and across model resolution. For instance, CRT for the Yellow Sea is overestimated by twofold in the 1°model but underestimated by 20-40% in the other versions. In the East Results are presented as the depth-than area-weighted, logarithmic mean of each cLMEs. Note that the coast-ocean boundary is extended ocean-ward by 1°for enhanced graphic presentation. (b-g) Comparisons of simulated CRTs between the 1°, 1/2°, and 1/4°model against the 1/8°model, respectively. Results are presented as the percentage difference relative to the 1/8°model. The red (blue) color is indicative of a positive (negative) bias in the coarser models relative to the 1/8°model. The gray color highlights the regions where model-model difference is insignificant (<25%). The RMSEs (root-mean-square errors) are calculated based on the log-transformed CRTs. Scatter size is scaled by the total volume of each cLME.
Bering Sea, however, simulated CRTs are all shorter than observations but successively improved with enhanced resolution (bias reduces from −57% to −25%).

Evaluation of Simulated Coastal Salinity Fields
A Taylor diagram shows that, globally, the 1/8°and 1/4°models have an overall high skill (r > 0.9, greater than the two coarser versions) with respect to capturing the observed coastal salinity deficits (Figures 2a and  2b), suggesting that coastal retention is better represented at high resolution. Figure 2c illustrates this with an example from the Northwest Atlantic shelf (cLME #7, 8, 9) where local bathymetry and circulation are complex ( Figure S5). The surface salinity field in this region is influenced by freshwater inputs from a large number of rivers (e.g., the St. Lawrence River), resulting in a shoreward decrease in surface salinity by approximately 5 units. The characterization of such a salinity gradient is poor in the 1°model but increasingly improved in models with enhanced resolution. A second example from the northern Gulf of Mexico is shown in Figure S6. Both examples show that freshwater plume structures are much better characterized at high resolution.

Sensitivity of Simulated CRT to Model Resolution
As the 1/8°model exhibits an overall higher skill in characterizing CRT, we used it as a baseline to assess the sensitivity of simulated CRT to model resolution (Figures 1b-1 g). When compared against the 1/8°model, the absolute percentage difference averaged across the 60 cLMEs decreases from 53% for the 1°model, to 37% for the 1/2°model, and to 26% for the 1/4°model. Specifically, coarser models are characterized by a negative CRT bias (indicative of too vigorous coast-ocean exchanges) in 39 of the 60 coastal systems (54% by area). This bias is particularly pronounced in narrow boundary currents. An example is the Oyashio Current in the western North Pacific, where CRT is underestimated by more than 50% in all three coarser models. However, this negative bias is not systematic across all cLMEs. For instance, a positive bias is seen in the Canadian Arctic Archipelagos and some semi-enclosed systems such as the Yellow Sea: Coarser models appear to create an unrealistically stagnant circulation in these regions.

The Geometric and Geostrophic Constraints on CRT
Overall, simulated CRT is longer on wide shelves at the northern high latitudes and lower in narrow boundary currents at low and the southern mid-latitudes (Figure 1). This pattern is broadly consistent with Sharples et al. (2017), which uses hydrodynamic principles to derive residence time of buoyancy plumes on the shelves as a function of latitude and shelf width. Specifically, Sharples et al. (2017) assume that the extent of a plume scales with the Coriolis parameter as a function of latitude (e.g., a plume at lower latitudes experiences a weaker Coriolis force hence would spread more quickly across The CRT estimator as a function of χ (the degree of geometric enclosure, χ = V/S) and ƒ (Coriolis frequency as a function of latitude) explains 73.2% (r = 0.856) of the system-to-system variability in simulated CRTs by the 1/8°model. Each filled circle represents one of the 60 cLMEs. Circle size is scaled by the total volume of each cLME. a shelf of a given width, resulting in a lower CRT). This estimator, however, does not consider the coastline shape (e.g., open versus semi-enclosed) or bathymetric variability that affects both the total coastal volume and the boundary area across which coastal and oceanic waters are exchanged. To account for these important factors, we tested the possibility of relating simulated CRT to the degree of geometric and geostrophic constraints. Specifically, we find that CRT (in days) for a given coastal domain can be expressed as: where a-c are constants (a = 1.198, b = 1.241*10 4 , and c = −8.339) and χ = V S is the ratio of the total volume (V in m 3 ) of a coastal system to its total open boundary area (S in m 2 ) across which coastal and oceanic waters are exchanged, calculated by summation of all grid cell's cross-sectional areas at the coast-ocean interface (Figure 3a). ƒ = 2Ωsinφ (in s −1 ) is the Coriolis frequency, with Ω the rotational rate of the Earth (Ω = 7.292*10 −5 rad s −1 ) and φ the median latitude (in rad) among all coastal grid cells within that system. Here the χ factor is indicative of the degree of geometric enclosure of the system, and ƒ relates to the Coriolis force that moves a cross-shelf flow toward the along-shelf direction, both positively correlated with CRT. This relationship explains 73.2% (r = 0.855) of the system-to-system variability in simulated CRT (Figure 3b). Additional regression analyses show that χ and ƒ factor alone explains 57.1% and 9.8% of the total variability, respectively ( Figure S10). This suggests that the degree of geometric enclosure plays a more important role in driving the variability of CRT than the geostrophic constraint.

Discussions and Conclusions
Coastal habitats are some of humanity's most important natural resources. It remains poorly understood how coastal ecosystems may respond to projected changes in the climate. This is in part constrained by our limited understanding of coastal retention timescales that may have significant implications on coastal ecosystem processes and material transport across the land-ocean interface (Monsen et al., 2002;Odebrecht et al., 2015). However, previous estimates are only available for a limited number of coastal sites and also strongly affected by the timing of sampling. We show that the distribution of CRT is highly heterogeneous on most shelves ( Figure S2) and that seasonal variability may exceed 40% (e.g., the southeast U.S. shelf; Figure S4). This suggests that residence time observed for a specific domain or site might not be representative for the entire shelf, and samples collected at specific times might not reflect the long-term mean properties in that system. Global syntheses on coastal residence time have been made only recently. Bourgeois et al. (2016) used a 1/2°g lobal ocean model to compute residence time for each specified coastal domain by dividing its volume by the integrated outflow at the coast-ocean boundary. Compared to our tracer-based method, this transportbased approach yields significantly shorter CRTs (p < 0.05; Text S3 and Figures S10 and S11). This difference may arise from how "coastal residence time" is defined. Specifically, our tracer-based method quantifies the time waters remain in any coastal domain (see details in section 2). Bourgeois et al. (2016), in contrast, assume that source waters to a specified domain have no previous "coastal history," and thus, the derived residence time only reflects how long those waters remain within that domain. While both perspectives have strengths and limitations, accounting for previous coastal history in our approach would yield longer residence times. This is illustrated for the Northwest Atlantic shelf, where alongshore currents with a long CRT serve as an important source to the region (Appendix S1 and Figure S8). Additionally, Bourgeois et al. (2016) heavily weight the rapid exchanges along the coast-ocean boundary. These exterior waters may exchange numerous times before the interior waters exchange once. Our approach, in contrast, averages residence times across the entire domain, providing a more distributed and thus generally longer characteristic CRT. The central Benguela Current provides a typical example, for which a strong cross-shelf CRT gradient is present ( Figure S7). Finally, it should be noted that, while developing a globally applicable approach is essential, alternative methods of deriving residence time exist and may be more appropriate for specific watersheds (e.g., Delhez et al., 2004 for river-dominated and semi-enclosed domains).
Compared to the 1/8°model, the three coarser models tend to produce shorter CRTs, indicative of more vigorous coast-ocean exchanges. For example, the negative CRT bias is pronounced in all boundary currents and also on most shelves in the mid-latitudes in the 1°model. This is in large part because coastal and open ocean grid cells are more directly connected in coarse models, prone to erroneous coastal flushing and an underestimated CRT (Text S2; Figures S7 and S8; Appendix S1). A few exceptions are noted in the Canadian Arctic Archipelagos and some semi-enclosed seas (i.e., Yellow Sea and the Gulf of Persian) where the coarser models tend to produce longer CRTs than the 1/8°model. In the case of the Canadian Arctic Archipelagos (cLME #18, 55, 66), coarsened model grids are unable to resolve the flows through small channels (usually represented by merely one or two grid cells), creating an unrealistically stagnant circulation locally and thus longer characteristic CRT for the system ( Figure S9). In 18 of the 60 coastal systems (38% by area), simulated CRTs are found consistent between the 1°and 1/8°models. These systems include coastal waters of the Southeast Asia, the Amazon shelf, the Mediterranean, the Somali Coastal Current, and several high-latitude shelves such as the Northern Bering-Chukchi Seas, the Laptev Sea, and the East Siberian Sea. This suggests that the representation of coastal circulation and retention is robust in these systems even at coarse resolution. Considering that most global ocean models (e.g., those participating in the Coupled Model Inter-comparison Project Version 6 or CMIP6) have a horizontal resolution of 1°to 1/4°, our comparison of simulated CRTs across model resolution may provide valuable insights into reliably synthesizing and interpreting shelf-scale information from multiple models at varying resolution.
We show that the degree of geometric enclosure (χ) and latitude of a coastal system together play a major role in driving CRT (equation ((1))). This provides a means for conveniently estimating CRTs for other interested coastal domains than those defined here (e.g., a segment of the southeast U.S. shelf as opposed to the entire cLME) and for evaluating the impact of coastline and bathymetric changes on CRT. As shown by the scattering in Figure 3b, the proposed CRT estimator is subject to biases in some coastal systems, for which other factors such as bathymetric heterogeneity may also play a critical role in determining CRT (e.g., two systems with a smooth vs. rough bottom may have a similar χ factor but a differing CRT due to small-scale circulation generated by bathymetric slopes). In addition, this estimator can only be related to the mean magnitude of the CRT and does not resolve its temporal variability as shown in Figure S4, for which changes in open ocean boundary circulation and freshwater inputs may also be important.
In summary, this study provides a novel, global perspective on coastal residence time (CRT) and its spatial and temporal variability across 60 coastal systems using a new tracer approach and a high-resolution (1/8°), state-of-the-art global ocean-ice model. Our model results agree well with previous observations made for a limited number of coastal systems. We also compared the CRTs simulated by the 1/8°model with three coarser models to discuss the sensitivity of simulated CRTs to model resolution. Our analysis shows that simulated CRTs have a negative bias (indicative of too vigorous coast-ocean exchanges) at coarser resolutions in 39 of the 60 coastal systems assessed (54% by area). This may represent an important mechanism underlying the poorly represented coastal ecosystem processes in coarse-grid models (e.g., the lack of coastal retention could prevent nutrients and other materials in source waters from impacting coastal phytoplankton biomass). Future incorporation of a biogeochemical component into our high-resolution models will lead us to a better understanding of carbon and nutrient transport from land to the ocean, and ultimately the role that coastal ocean plays in global biogeochemical cycles.