Overview of (Sub)mesoscale Ocean Dynamics for the NAAMES Field Program

TheNorth Atlantic Aerosols andMarine Ecosystems Study (NAAMES) is an interdisciplinary effort to characterize plankton ecosystem properties through the annual cycle and determine how remote marine aerosols and boundary layer clouds are influenced by marine ecosystems, especially by phytoplankton (Behrenfeld et al., 2018). This study was carried out in the nortwestern Atlantic, a region characterized by intense mesoscale and submesoscale variability that has been observed to affect phytoplankton abundance, phenology and community composition (McGillicuddy, 2016).


INTRODUCTION
The North Atlantic Aerosols and Marine Ecosystems Study (NAAMES) is an interdisciplinary effort to characterize plankton ecosystem properties through the annual cycle and determine how remote marine aerosols and boundary layer clouds are influenced by marine ecosystems, especially by phytoplankton (Behrenfeld et al., 2018). This study was carried out in the nortwestern Atlantic, a region characterized by intense mesoscale and submesoscale variability that has been observed to affect phytoplankton abundance, phenology and community composition (McGillicuddy, 2016).

(Sub)mesoscale Ocean Variability
The oceanic (sub)mesoscale (defined here by spatial scales of O(1 − 100 km) and temporal scales from days to months) encompasses dynamical processes and features such as mesoscale meanders, eddies, and fronts, as well as submesoscale fronts and vortices. These features affect phytoplankton by structuring its distribution, its access to resources (i.e., nutrients, light), the competition to exploit these resources, and its encounter rates with grazers and viruses (i.e., top-down control). In this section we summarize the different mechanisms by which (sub)mesoscale features can influence phytoplankton and thus play a role in the creation of biogenic aerosols.

Mesoscale Eddies and Meanders O(10 − 100 km)
Bottom-up effects of mesoscale eddies and meanders on phytoplankton vary regionally throughout the global ocean (Gaube et al., 2014). In eddies, vertical fluxes are driven by the displacement of isopycnals either by internal dynamics (e.g., eddy-pumping or eddy/eddy interaction) or as the result of a surface forcing (e.g., eddy-induced Ekman pumping and eddy-mediated changes to mixing). The intensification of cyclonic eddies results in the upward displacement of isopycnals and upwelling of nutrients into the euphotic zone. Conversely, a downward displacement of isopycnals results in the downwelling of nutrients and phytoplankton during the intensification of anticyclonic eddies. Ekman pumping in the center persists for the entire lifetimes of cyclones and anticyclones and is forced by the interaction of rotating surface currents with the wind. Eddy-induced Ekman pumping results in upwelling in the cores of anticyclones (Dewar and Flierl, 1987) and downwelling in the cores of cyclones (McGillicuddy et al., 2007;Gaube et al., 2015). Furthermore, anticyclones and cyclones are generally characterized by deeper and shallower mixed layers compared to the surrounding waters (Dewar, 1986;Dufois et al., 2014;Hausmann et al., 2017;Gaube et al., 2018). Deeper mixed layers in anticyclones can result in nutrient fluxes that would not happen in cyclones in the same region (Gaube et al., 2013;Dufois et al., 2014). In addition, these differences in mixed layer depth are likely to dilute or concentrate phytoplankton and their grazer and therefore decouple or couple these populations, modulating trophic interactions.
The horizontal advection of nutrients and plankton can be summarized in two distinct mechanisms: (1) stirring, which occurs primarily around the peripheries of eddies; and (2) trapping and subsequent transport in the interiors of eddies of both phytoplankton and their grazers (Isla et al., 2004;Gaube et al., 2014). In nonlinear eddies, rotational velocities exceed the eddy propagation rate, trapping fluid in their interiors and creating shear zones along the peripheries of eddies (McWilliams and Flierl, 1979;Flierl, 1981). As a result, ecosystems trapped in eddies during formation act to "prime" the eddy interior toward either elevated or suppressed phytoplankton concentration and can potentially contain planktonic communities that are distinct from their surroundings (Bracco et al., 2000;d'Ovidio et al., 2010).
The interplay of this mosaic of mechanisms by which eddies modulate phytoplankton communities is complex and varies according to the region of the origin of eddies, the areas into which they propagate, and their life stage (Gaube et al., 2014;Frenger et al., 2018). In the North Atlantic, the trapping and subsequent transport of phytoplankton by eddies is hypothesized to be the dominant bottom-up mechanism by which eddies influence phytoplankton (Gaube et al., 2014). However, there is observational support for the potential impacts of eddy-mediated vertical nutrient fluxes on phytoplankton within eddy interiors (Franks et al., 1986;Gaube and McGillicuddy, 2017).
Importantly, eddies do not only affect the quantity and quality of phytoplankton, but they also influence the timing and seasonal progression of bulk phytoplankton properties. In the North Atlantic, the world's largest spring phytoplankton bloom results in peak chlorophyll in the late spring through early summer, but the precise timing of the maximum has been shown to be influenced by the presence of mixed-layer eddies (Mahadevan, 2016). In oligotrophic regions, eddies have also been shown to influence the timing of peak chlorophyll. For example, in the South Indian Ocean, an eddy-centric analysis constructed from thousands of long-lived eddies showed that the seasonal maximum in chlorophyll occurred nearly a half a month earlier in anticyclones when compared to the surrounding region (Gaube et al., 2014). The influence of eddies on the seasonal timing of the chlorophyll maximum appears to be the result of eddy influences on mixed-layer depth and the mixing of nutrients into the euphotic zone (Gaube, 2012;Dufois et al., 2014;Gaube et al., 2014;Mahadevan, 2016).

Submesoscale Eddies and Fronts O(1 − 10 km)
Submesoscale buoyancy gradients, or fronts, arise in the surface ocean due to instability (often resulting from the strain and stirring generated by mesoscale eddies) or surface forcing such as wind stress and fluxes of heat and freshwater into and out of the ocean surface. Submesoscale fronts are associated with strong vertical velocities in the mixed layer resulting in the transports of nutrients and phytoplankton between the surface and the thermocline. This has implications for the structure of phytoplankton communities in terms of modulation of the ambient light or nutrient field (Lévy et al., 2012(Lévy et al., , 2018Mahadevan, 2016). Furthermore, submesoscale fronts are regions of convergence/divergence and dilution induced by (sub)mesoscale currents. This can have important implications for the encounter rate between phytoplankton and their grazers and can result in accumulation of phytoplankton biomass (Calil and Richards, 2010;Lehahn et al., 2017). Linking these physical processes to their biological response is necessary to explain the spatial heterogeneity in phytoplankton biomass and community structure, and hence the creation of biogenic aerosols. While altimetry cannot resolve explicitly the submesoscale, a Lagrangian re-analysis of these field can reveal submesoscale features induced by the stirring caused by mesoscale eddies (Hernández-Carrasco et al., 2011).

Using Satellite Altimetry to Aid in the Interpretation of Airborne and Ship-Based Observations
To estimate the location and properties of mesoscale features we used the sea level anomaly (SLA), mean dynamic topography (MDT -defined as the 20-year average of sea surface height measured by a constellation of satellite altimeters) and geostrophic currents distributed by the Copernicus Marine Environment Monitoring Service (CMEMS, http://marine. copernicus.eu). The altimetry-derived geostrophic velocities used to compute the Lagrangian diagnostics detailed in section 4 are estimated from the absolute dynamic topography (ADT) which is defined as the sum of SLA and MDT. Mesoscale eddies were identified and characterized from high-pass filtered SLA fields. The high-pass filtered SLA observations are defined as the difference between the daily SLA fields and the smoothed SLA that captures the effects of seasonal heating and cooling following the procedure laid out in Chelton et al. (2011b). These high-pass filtered SLA fields are simply referred to as SLA in this manuscript and are structured by mesoscale features such as eddies and meanders. The MDT product was used to define a subregional classification for the different locations explored during the field program (section 2). For the 2015 to 2017 expeditions, we used the Delayed Time products (DT, former UPD) and for the 2018 expedition we used the Near Real Time product (NRT). Both products are distributed on a daily basis on 1/4 • grids.
Areas within the study region that are expected to be influenced by mesoscale eddies and meanders can be estimated by quantifying the portion of time a given location is inside a mesoscale eddy or current meander (defined collectively as coherent mesoscale structures, or CMS). This measure of eddy coverage varies nearly an order of magnitude over the greater NAAMES regions (Figure 1), with nearly 80% coverage in the proximity of the Gulf Stream to less than 10% in the location of the northernmost NAAMES stations. In light of the spatial heterogeneity of eddy coverage, amplitude, and size, along with the variety of mechanism by which eddies impact biological and biogeochemical processes, we provide a brief analysis of the mesoscale eddy field at each NAAMES stations in the Supplementary Materials of this article.
The data presented here are publicly available in NetCDF format as part of the SeaWiFS Bio-optical Archive and Storage System (SeaBASS) maintained by the NASA Ocean Biology  Processing Group and can be accessed at the URL https://seabass. gsfc.nasa.gov/archive/UWASH/gaube/NAAMES/documents.

A SUBREGIONAL CLASSIFICATION SCHEME FOR THE NORTH ATLANTIC
The NAAMES expeditions covered a large region of the North Atlantic with significant spatial variability in physical and biological properties (Behrenfeld et al., 2018). To first order, much of this variability scales with latitude. A simple geographic binning of the NAAMES observations, however, would not correctly capture the spatial variability in the mesoscale eddy field (Figure 1). We propose a subregional classification scheme that is based on the MDT (Figure 2). We have chosen four subregions by visually separating the MDT map into the most obvious different areas: (1) subarctic, defined as regions with MDT ≤ −51 cm; (2) temperate, defined as the area with MDT in the range −51 cm < MDT ≤ −10 cm; (3) subtropical, defined as the region with MDT in the range −10 cm < MDT ≤ 30 cm; (4) and the Sargasso Sea and Gulf Stream (MDT > 30 cm). These regions each represent distinct physical provinces and can be characterized by a general southward gradient in temperature, resulting in a homologous gradient in MDT.

OVERVIEW OF COHERENT MESOSCALE STRUCTURES IN THE NAAMES REGION
The northwestern Atlantic is a region that encompasses large amplitude CMSs, with its largest features falling in the upper 95th percentile of amplitude globally (Chelton et al., 2011a). The largest amplitude CMSs are found in the Gulf Stream region where average amplitude can exceed 40 cm ( Figure S1). These features are generated as mesoscale meanders of the Gulf Stream becoming unstable and pinching off to become eddies. In the region of focus for NAAMES (see black rectangle in Figure S1), there is a general northward gradient in decreasing CMS amplitude, with a region of moderate amplitude CMS (≈ 20 cm) flowing from the east near the middle of the domain. Eddies in the northern portion of the NAAMES region (including the Temperate and Subpolar regions, Figure 2) likely form in the open ocean as a result of baroclinic instabilities.
The distribution of CMS amplitudes in the NAAMES region is not symmetric with respect to polarity; a larger proportion of large amplitude CMSs are cyclones, when compared to anticyclones ( Figure S2a). Conversely, a greater number of large radius eddies are anticyclones ( Figure S2b). This combines to result in higher intensity CMS, defined here as the ratio of amplitude to radius, that are preferentially cyclonic (Figure S2c).

LAGRANGIAN RE-ANALYSES OF ALTIMETRY
Lagrangian diagnostics are generally defined as quantities calculated in the frame of reference of a water parcel that is followed through time (Kundu et al., 2008;van Sebille et al., 2018). In the context of NAAMES, Lagrangian diagnostics can be helpful to identify frontal regions, eddies cores and peripheries and to estimate where the sampled water parcels are coming from. For this specific study we used the LAMTA Lagrangian scheme (available on bitbucket.org/f_nencio/spasso/overview) to advect water parcels using altimetry-derived current velocities (d' Ovidio et al., 2015). Since estimating horizontal velocities from altimetry requires making the assumption of geostrophic equilibrium, the velocity field and the derived Lagrangian diagnostics are expected to describe the ocean geostrophic layer and neglect the top surface layer where wind driven circulation is dominant (Ekman layer).

Caveats
1. The Lagrangian diagnostics are computed on a highresolution grid (1 km grid spacing), yet the calculation of the velocities assumes a geostrophic balance and therefore do not include any "unbalanced" motions (such as wind-driven currents). Therefore this method only identifies submesoscale features that are generated by the currents at the mesoscale and larger (Keating et al., 2012). 2. All the Lagrangian diagnostics presented in this study are calculated from horizontal velocities. This means that if a tracer has been transported to the surface by any kind of vertical movement, the diagnostics may not correctly identify its origin. 3. Even at a higher resolution, the geostrophic velocity field still defines a chaotic system. As a consequence it is impossible to calculate the "true" trajectories of specified parcels. Errors in the definition of the initial location of a water parcel and in the velocity field increase with increasing time of advection. Therefore the calculated trajectories can only be uniformly close to a true trajectory (Özgökmen et al., 2000;Prants et al., 2018).

Finite Size Lyapunov Exponents-an Index of Frontal Activity [Units = d −1 ]
Finite Size Lyapunov Exponents (FSLE) are defined as the rate of confluence of water parcels initially far apart (Aurell et al., 1997;d'Ovidio et al., 2004). High values of FSLE refer to regions where water parcels that were initially far have converged rapidly. FSLE are defined as: (lon, lat, t) refer to the location in time and space of the region where the FSLE is computed, τ is the time taken for water parcels at an initial distance of δ (in this study = 0.3 • ) to reach a final separation of δ 0 (in this study = 0.01 • ). In this analysis τ can reach a maximum value of 60 days. Evaluation of the robustness of FSLE can be found in Cotté et al. (2011) andHernández-Carrasco et al. (2011).

Eddy Retention Parameter-How Long Has a Water Parcel Been in an Eddy [Units = d]
The eddy retention parameter (RP) quantifies for how long a water parcel has been recirculating within the core of an eddy (d' Ovidio et al., 2013) defined as a region of negative values of the Okubo-Weiss parameter (OW) (Okubo, 1970;Weiss, 1991) [units = d −2 ]. The OW is calculated as: OW(t, lon, lat) = s 2 (t, lon, lat) − ω 2 (t, lon, lat); (2) (3) where x and y refer to longitude and latitude and v x and v y are the zonal and meridional components of the velocity field. The OW parameter represents the difference between the vorticity of the velocity field and the strain. It is expected to be negative in regions where the rotational component of the velocity field dominates over the strain (e.g., within eddy cores) and positive in regions where strain dominates over vorticity (e.g., at the peripheries of eddies). Given a water parcel and its back-trajectory, the RP parameter accounts for how many days prior to the sampling date a given water parcel has been located in a region with negative OW (i.e., within an eddy core). In this study the maximum time used for the backward advection of water parcels is 30 days, thus the most retentive features have RP = 30 days.

Origin of Water Parcels-Where Do Tracers Come From? [units = • ]
These diagnostics are a measure of the latitude and longitude of water parcels 15 days before the sampling date. This advection time was chosen as a compromise between timescales of interest for phytoplankton ecology, computational costs, and the limitations due to the accumulation of error mentioned section 4.1.

Water Origin Mask Fields [Units = d]
In order to visualize the origin of the water parcels sampled by the R/V Atlantis, we developed a masks product that indicates the locations of the water parcels within a given radius from the ship position as a function of time prior of being sampled. Such locations are labeled by the number of days prior to when the ship interacted with that specific water parcel (up to 15 days). All remaining locations are flagged. The mask were calculated every 6 hours for the NAAMES region with a spatial resolution of 0.05 • .
The masks were computed as follows: 1. The location of R/V Atlantis at a given time and date was used to define the center of two disks having 20 and 40 km radii (see Figure S3a). 2. All the water parcels included within the disks were advected backward in time for 15 days using the altimetry-derived velocities described in section 4 (see Figure S3b which shows the trajectories of individual water parcels at a given time using the 40 km disk). 3. Using the trajectories from (2) we created a mask file that indicates how many days prior each 0.05 • grid point interacted with the water parcel defined by the disks described in (1). This results in individual masks for both the 20 km ( Figure S4a) and 40 km (Figure S4b) disks. Table 1 and Figure S5 summarize the stations sampled during the four NAAMES expeditions and their properties in respect to the mesoscale eddy field. Out of the total of 33 stations sampled in the course of the four expeditions, one was located within a mode-water eddy, seven within anticyclonic features and Red lines refer to anticyclonic eddies, blue lines to cyclonic eddies and the yellow line refers to the mode-water eddy sampled during NAAMES3. Solid red and blue refer to eddy cores and light red and blue refer to eddy peripheries. (*) The dates refer to the arrival on station. (**) Retentive surfaces correspond to the largest closed surface around an eddy presenting a retention higher than 4 days.

SUMMARY OF EDDY FEATURES SAMPLED DURING THE NAAMES PROGRAM
seven within cyclones. The seasonal and regional distribution of the sampled eddies is not uniform: the choice of the NAAMES stations across four expeditions was a compromise between covering a large latitudinal gradient, sampling eddies having different properties, tracking bio-argo floats to create a multi-year time series of physical bio-optical parameters and mitigating the operational challenges of working in the North Atlantic. During NAAMES1, corresponding to the winter transition phase of the seasonal phytoplankton cycle (Behrenfeld et al., 2018), the peripheries of three anticyclonic eddies and one cyclonic eddy as well as the core of a cyclonic eddy were sampled. In the following spring, during the climax transition phase, two anticyclones' peripheries, one anticyclone's core, one cyclone's periphery and one cyclone's core were sampled.
During NAAMES3, corresponding to the declining phase, we sampled an anticyclone's core, a cyclone's core, a mode-water eddy's core and a cyclone's periphery. During NAAMES4, corresponding to the accumulation phase only the periphery of a cyclonic eddy was sampled. Across different seasons, we sampled a cyclone and an anticyclone in the subpolar subregion, a cyclone and two anticyclones in the temperate subregion, and an anticyclone in the Gulf Stream and Sargasso Sea subregion. The most extensively sampled subregion was the subtropical one, with stations within five cyclones, three anticyclones and a mode-water eddy were sampled as well as six out of eddy stations. The eddy properties of each station mentioned in Table 1 are detailed in the Supplementary Material (Figures S6, S7 and following).

AUTHOR CONTRIBUTIONS
AD and PG contributed equally to the compilation of the data, their analysis and interpretation and the writing of the manuscript.

FUNDING
This study has been conducted using E.U. Copernicus Marine Service Information (CMEMS). AD is grateful for the support of the Applied Physics Laboratory Science and Engineering Enrichment Development (SEED) fellowship. This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 749591. PG is grateful for the support of the NASA Physical Oceanography Program (NASA grant NNX16AH59G) and the NSF award OCE-1558809.