Seasonality of submesoscale flows in the ocean surface boundary layer

A signature of submesoscale flows in the upper ocean is skewness in the distribution of relative vorticity. Expected to result for high Rossby number flows, such skewness has implications for mixing, dissipation, and stratification within the upper ocean. An array of moorings deployed in the Northeast Atlantic for 1 year as part of the experiment of the Ocean Surface Mixing, Ocean Submesoscale Interaction Study (OSMOSIS) reveals that relative vorticity is positively skewed during winter even though the scale of the Rossby number is less than 0.5. Furthermore, this skewness is reduced to zero during spring and autumn. There is also evidence of modest seasonal variations in the gradient Rossby number. The proposed mechanism by which relative vorticity is skewed is that the ratio of lateral to vertical buoyancy gradients, as summarized by the inverse gradient Richardson number, restricts its range during winter but less so at other times of the year. These results support recent observations and model simulations suggesting that the upper ocean is host to a seasonal cycle in submesoscale turbulence.


Introduction
The ocean surface boundary layer (OSBL) is an important component of Earth's climate system. It serves as the medium through which heat, momentum, and gases (e.g., carbon dioxide) are exchanged between the ocean and atmosphere through a combination of physical, chemical, and biological processes. Prediction of climate on seasonal, decadal, and longer time scales therefore requires an accurate knowledge of the dynamics within this layer.
Wind, waves, precipitation, and heat loss at the ocean surface have long been known to have a significant effect on vertical stratification within the OSBL and have successfully been modeled as one-dimensional processes [see Large et al., 1994, and references therein]. More recently, however, a host of three-dimensional processes have been identified as influencing both stratification and dissipation in the upper ocean. These include (1) effects of Coriolis-Stokes forcing and Langmuir turbulence, (2) circulation due to frontogenesis/ frontolysis, and (3) forced and unforced instabilities [Hoskins, 1982;Thomas, 2005;Polton et al., 2005;Polton and Belcher, 2007;Boccaletti et al., 2007;Thomas et al., 2008;Fox-Kemper et al., 2008a;Taylor and Ferrari, 2009;Belcher et al., 2012;D'Asaro et al., 2011;Hamlington et al., 2014;Brüggemann and Eden, 2015]. The dynamical regime in which these processes reside is the oceanic submesoscale.

Skewness in the Distribution of Relative Vorticity
One of the hallmarks of submesoscale flows is skewness in the distribution of relative vorticity, . This skewness is more common for submesoscale flows compared with motions at the mesoscale (30-200 km) since the Rossby number is considerably smaller at the mesoscale. Both Rudnick [2001] and Shcherbina et al. [2013] document positive skewness in winter from observed distributions of estimated from shipboard velocity measurements, and similar skewness has been reported by Roullet and Klein [2010] and Qiu et al. [2014] in high-resolution ocean models. The simplest dynamical explanation for this observed asymmetry is that small perturbations grow exponentially when potential vorticity (PV) takes on sign opposite the Coriolis parameter, is the Ertel PV [Ertel, 1942;Hoskins, 1974]. Here (f + ) is the vertical component of absolute vorticity, h is the horizontal vorticity vector, N 2 is defined above, and ∇ h b is the horizontal buoyancy gradient. Because relative vorticity is one factor determining q and since instabilities tend to drive fq toward a positive state, these instabilities restrict the range of observed values of . In particular, extreme negative values of are rare in the Northern Hemisphere, and positive values of are rare in the Southern Hemisphere. This results in a positively skewed probability density function (PDF) of Ro = ∕f . Note that increased variance in Ro resulting from vortex stretching will broaden the distribution of Ro and help reveal this skewness [Capet et al., 2008].

Focus of This Study
There are few open ocean measurements of submesoscale Ro. Those of Rudnick [2001] are biased low owing to an approximation of from along-track gradients alone, and Shcherbina et al. [2013] document Ro in a region influenced by the mean flow of the Gulf Stream. Furthermore, in most cases, existing observations are limited to wintertime conditions in which, based on modeling [Mensa et al., 2013;Qiu et al., 2014;Sasaki et al., 2014;Brannigan et al., 2015] and observational [Ostrovskii, 1995;Callies et al., 2015] studies, we expect submesoscale flows to be most active.
In this study, we extend the observational efforts of Rudnick [2001] and Shcherbina et al. [2013] by estimating Ro = ∕f from measurements spanning an annual cycle. We conduct this analysis on observations from the Northeast Atlantic, a region characterized by weak mean kinetic energy and low-to-moderate mesoscale eddy kinetic energy ( Figure S1 in the supporting information). We demonstrate that relative vorticity within the OSBL is positively skewed during winter even though a characteristic scale for |Ro| is less than 0.5 at all times. Additionally, this skewness appears to vary seasonally, a result that has implications for stability. In summary, we (i) present rare observational evidence of submesoscale flows in the open ocean, (ii) provide a scale for |Ro|, and (iii) highlight seasonal changes in turbulence parameters with implications for stability.  Figure 1a) [Allen et al., 2013]. Measurements include temperature, salinity, and horizontal velocity from moorings (50-500 m) and hydrographic observations from two ocean gliders (0-1000 m). The present study focuses on mooring observations. Two clusters of moorings were oriented around a central mooring. One cluster consisted of four moorings spaced approximately 1.5 km from this center mooring, while the other cluster consisted of four moorings spaced approximately 10 km from this same mooring ( Figure 1b). We refer to regions spanned by these moorings as inner and outer mooring domains, respectively, and focus attention on inner moorings.

Observations Within the OSBL
MicroCAT conductivity-temperature-depth (CTD) sensors and acoustic current meters (ACMs) were spaced approximately 30-100 m in the vertical, with increased numbers on the center mooring in order to resolve stratification and shear ( Figure S2). Sampling intervals of CTDs and ACMs were 5 and 10 min, respectively. Velocity, temperature, and salinity at each mooring were interpolated to common depths and times, and smoothing was applied in the vertical and time (30 m and 2 h) in order to minimize noise while retaining time and spatial scales of interest. In an effort to prevent damage by severe weather throughout the year, OSMOSIS moorings resided below the ocean surface. This precluded use of a Global Positioning System (GPS), and owing to their considerable cost, acoustic positioning sensors were not used. For these reasons, knowledge of mooring position is somewhat limited. In the work that follows, we compute vorticity from neighboring moorings and, therefore, require accurate position information. To quantify the error introduced by position uncertainty, we modeled the distance between moorings as the sum of fixed distances and a stochastic parameter that represents perturbations from these constant distances. Fixed distances were estimated from anchor positions and perturbations were modeled as random variables with zero mean and standard deviations determined from the time integration of horizontal currents (Appendix A). A Monte Carlo approach was then used to quantify the range of uncertainty. Such an approach provides an upper bound on errors due to mooring motion since it effectively assumes instruments drift as unconstrained floats.

Vorticity From the Inner Moorings
We divided the inner mooring domain into four triangular regions in order to explore submesoscale characteristics described above. Relative vorticity was estimated from circulation following Stokes' theorem [Batchelor, 2000]: Here A is the area enclosed by the path, u h is the horizontal fluid velocity, and dl is a differential pointing in the direction of the closed path. For the present case, our closed path is the area formed by a single triangle (cf. Figure 1b). We estimated the horizontal velocity, u h , as the average of velocities between two neighboring moorings and |dl| as the distance between these moorings. Alternatively, the vertical component of relative vorticity can be estimated from velocity gradients, yielding numerically equivalent values [Bryden and Fofonoff , 1977].
From the above expression, Ro = ∕f was computed (f =1.094 × 10 −4 s −1 ), and the standard deviation and skewness of Ro were estimated as functions of depth and season. The choice of depth and season are discussed below. We combined vorticity estimates from all four inner triangles in order to obtain robust statistics. We also computed seasonal profiles of stratification and gradient Richardson number from buoyancy and shear at the center mooring. Note that as Ri is approximately lognormal, averaging must be done in logged coordinates (before transforming it back to a characteristic value). Also, because the mean of Ro is approximately zero, the standard deviation of Ro provides a characteristic scale for the Rossby number.
We have included estimates of vertical stratification from ocean gliders [Thompson et al., 2016]

Defining Depths and Seasons
Mooring measurements made within the OSBL and where we expect submesoscales to be most active do not extend for the entire observational record. This is most easily discerned from two-dimensional histograms of observation number as a function of time and normalized depth, |z|∕h (not shown). Here |z| is the observation depth and h is the mixed layer (ML) depth, defined as the depth at which potential density exceeds that observed at 20 m by 0.02 kg m −3 [de Boyer Montégut et al., 2004]. We used these histograms to define the following two seasons: winter (20 December to 23 April) and autumn/spring (hereafter, nonwinter; 11 October to 19 December and 24 April to 3 June). Additionally, to facilitate comparison across seasons, we have computed statistics in |z|∕h coordinates. Here h has been smoothed using a 12 day window to eliminate intermittent changes in ML depth. With a few exceptions (e.g., early February and May), h tracks the depth of the pycnocline well (Figure 1c). We argue that this depth most appropriately characterizes submesoscale vortices since baroclinic instability (BCI) of the surface ocean is inherently tied to a local but nonintermittent index of ML depth [Stone, 1966;Boccaletti et al., 2007;Fox-Kemper et al., 2008b;Grooms, 2015].

Results
Year-long records of stratification, N∕f , from ocean gliders and normalized relative vorticity, ∕f , within the inner mooring domain are shown in Figures 1c and 1d. Overlaid on these are ML depth and contours of constant |z|∕h. Note that Ro has values on the order of 1 but is generally smaller. Additionally, a number of short-duration (i.e., < 3 days) anomalies propagate through the inner mooring domain over the observational period, suggesting a sufficient number of samples for statistical estimates.
In Figure 2, we display statistics of stratification, vorticity, and gradient Richardson number. We also highlight the depth of peak stratification, |z|∕h = 1.3, and refer to locations above this depth as being within the OSBL.
Representative stratification profiles depict a strong seasonal cycle (Figure 2a). Nonwinter values are as high as N∕f = 40, while winter values fall below N∕f = 15. This pattern is reproduced by glider measurements. Higher nonwinter magnitudes are present within the OSBL owing to the glider's ability to measure closer to the ocean surface.
Standard deviations of ∕f also reveal seasonal differences (Figure 2b). Differences are pronounced within the OSBL, but there is still evidence of seasonality immediately below the OSBL. Standard deviations of ∕f within the OSBL are approximately 0.4 during winter and 0.3 during nonwinter. These serve as scales for the Rossby number. Also noteworthy, the standard deviation of ∕f decreases with increasing depth during both periods, supporting concepts of surface-intensified flow.
Similarly, there are seasonal differences observed in the skewness of ∕f (Figure 2c). Skewness within the OSBL exceeds 1.0 during winter but is virtually nonexistent during nonwinter. The existence of positive skewness at depth below the OSBL during nonwinter can be traced to anomalously strong cyclonic flow below the OSBL during late November and early December (cf. Figure 1d).
Finally, the logged-gradient Richardson number also varies seasonally (Figure 2d), with winter magnitudes close to 2 and nonwinter magnitudes of approximately 6. Such seasonality has also been documented in a recent study involving ocean gliders [Thompson et al., 2016].

Seasonality in Vorticity Skewness
The distribution of Ro = ∕f is positively skewed within the OSBL during winter despite that a scale for the Rossby number is 0.5. What gives rise to vorticity skewness if, strictly speaking, |Ro| is not O(1)? To address this question, let us revisit the stability criterion. First, in the absence of lateral buoyancy gradients and assuming positive stratification, the condition for stable flow becomes, fq > 0 →f 2 (1 + Ro)N 2 > 0, or Ro > −1. This is the well-known bound on Ro for centrifugally stable flow. Second, assuming thermal wind balance and w ≪ u, v, one can rewrite equation (1) as In this case, the condition for stable flow simplifies to Ro > −1 + Ri −1 , where Ri = f 2 N 2 ∕|∇ h b| 2 is the balanced gradient Richardson number. Instability occurring for values less than −1 + Ri −1 is referred to as symmetric instability [Hoskins, 1974;Colin de Verdière, 2012]. This expression states that the PDF of is influenced not only by its value relative to −f but also by the ratio of lateral and vertical buoyancy gradients. This could explain why we observe vorticity skewness for only moderate |Ro|. Ignoring for the moment small seasonal changes in Ro, it also provides a plausible explanation for the observed seasonality in vorticity skewness, as illustrated in Figure 3.

Discussion
The ocean surface is subject to atmospheric forcing in the form of winds, precipitation, and heat. Such forcing, particularly strong during winter, alters PV and can drive it toward zero [Thomas and Lee, 2005;Thomas, 2005]. In this case, q is not conserved and a host of instabilities can develop [Stone, 1966;Hoskins, 1974;Haine and Marshall, 1998;Taylor and Ferrari, 2009;Thomas et al., 2013;Lazar et al., 2013]. Note that BCI [Boccaletti et al., 2007] is not included here since it can occur for positive fq [Stone, 1966].
The results presented above depict an active OSBL consistent with this interpretation. Decreased stratification in winter permits enhanced vertical motion. Through a combination of vortex stretching due to frontogenesis [Hoskins, 1982], stretching due to interfacial waves at the base of the OSBL and vortex generation by baroclinic and shear instabilities [Stone, 1966;Munk et al., 2000;Eldevik and Dysthe, 2002;Boccaletti et al., 2007;Grooms, 2015], the standard deviation of Ro = ∕f increases. Finally, as |Ro| increases and Ri decreases, the lower bound associated with symmetric instability emerges.
Increased variance in Ro due to vortex stretching has the effect of revealing vorticity skewness even in the absence of instabilities [Capet et al., 2008]. Thus, whether or not the marginally stable state (i.e., fq = 0) is crossed, leading to the suite of instabilities referred to above, and to what degree wind forcing and buoyancy loss at the surface are relevant factors [Thomas and Lee, 2005;Thomas, 2005;Taylor and Ferrari, 2010;Thomas and Taylor, 2014;Hamlington et al., 2014;Whitt and Thomas, 2005] are questions that we leave for a future study.

Limitations and Sources of Error
Results of this study depend upon observations within a single year (September 2012 to September 2013). Furthermore, owing to the minimum observable depth of moorings, the number of measurements made throughout the year is restricted to those in late Autumn through early Spring; that is, we are unable to comment on summertime statistics. While intermittency presents a challenge, we can nevertheless estimate confidence levels on statistics. To do this, the effective sample number, n, was estimated as the total duration of time a given normalized depth was observed divided by t d = 3 days. (The decorrelation time scale, t d , of u and v is 1 day, but we have opted for a more conservative value.) Standard errors of each statistic were then computed and used to determine the significance of results. In all cases, the results given above remain unchanged ( Figure S3).
Another limitation pertains to the mooring array. Owing to the fixed spacing of moorings, we are unable to distinguish between (a) seasonally varying amplitudes and (b) seasonally varying scales of submesoscale phenomena. In particular, during winter, the ML deformation radius is approximately d = (N∕f )h = 4 km, while during summer months, it can be closer to d = 1-2 km. In fact, this is a problem inherent in all studies thus far [Mensa et al., 2013;Qiu et al., 2014;Sasaki et al., 2014;Brannigan et al., 2015;Callies et al., 2015] as their horizontal spacings only marginally resolve the seasonally varying deformation radius.
Uncertainty has been introduced into this study as a result of mooring motion. Such motion gives rise to two types of errors. The first is additional variance inserted into ∕f . The second is that observed velocities will be smaller than true velocities if moorings are advected with the flow. Consequently, vorticity (cf. equation (2)) will be biased low and standard deviations of ∕f reported here might be biased low. Note, however, that these two types of errors tend to compensate one another.

Conclusions
There is increasing evidence from modeling studies [Mensa et al., 2013;Qiu et al., 2014;Sasaki et al., 2014;Brannigan et al., 2015], and observational studies [Ostrovskii, 1995;Callies et al., 2015;Thompson et al., 2016] that submesoscale flows undergo a seasonal cycle. The observations presented here support this concept by depicting seasonal changes in gradient Rossby and Richardson numbers within the OSBL: (i) variances of Ro are enhanced during winter months over nonwinter months, (ii) Ro is positively skewed during winter months while approximately zero at other times of the year, and (iii) we observe seasonal changes in Ri, consistent with seasonal changes in vorticity skewness. These observations have potential implications for mixing, stratification, and dissipation within the OSBL with application for modelers seeking to represent submesoscale processes in climate-scale ocean models. Note, for example, that symmetric instability has yet to be parameterized in these models.
An important finding of this study is that submesoscale |Ro| is approximately 0.3-0.45 in the Northeast Atlantic. These values are smaller than those observed near the Gulf Stream (e.g., ∼0.7-1.0) [Shcherbina et al., 2013] and suggest a different balance of terms in the momentum equations. This may find wider application, for example, in other open ocean regions where both mean and eddy kinetic energies are smaller.

Appendix A: Mooring Motion
As mentioned in section 2.1, we have limited knowledge of the relative positions of moorings in time. Here we quantify the error produced by such uncertainty using a Monte Carlo simulation.

A1. Mooring Design
Considerable buoyancy was placed on mooring cables between 0 and 500 m depth so that the moorings stood as close to vertical as possible within the OSBL. This created an effective pivot point within the catenary at 600 m as evidenced by mooring design simulations (P. Provost, personal communication, 2014). In particular, the depth-integrated flow determines the average positions of moorings, but surface-intensified flows appear to perturb the positions of moorings with respect to one another.

A2. Coherence of Motion at the Inner Moorings
Under the assumption that surface currents determine the relative position of moorings, we estimated the coherence and phase of u h between the center and each of the inner moorings ( Figures S4a and S4b). Horizontal velocities in the upper ocean are coherent for time scales greater than 12 h (i.e., mean squared coherence, | | 2 ≥ 0.6) with motions less than 6 h being considerably incoherent (i.e., | | 2 <0.2). This is also reflected in the phase, which becomes noisy for periods less than 6 h but is well-defined for longer periods. While in part consisting of fluctuations from wind and waves, these shorter-period motions contain the signal of interest. Low-pass filtering the signal is therefore imprudent, and we instead model relative mooring motion as a stochastic process.

A3. A Stochastic Model for Mooring Motion
The relative position of moorings is expressed as the sum of a fixed distance and a perturbation distance. Fixed distances were determined from estimated anchor positions, while perturbation distances, = ( x , y ), were modeled as Gaussian white noise processes with zero mean and nonzero variance, 2 = ( 2 x , 2 y ), estimated from the time integration of differential currents. The distance between center and neighboring moorings are assumed uncorrelated, and motion is constant with depth at a given time.

10.1002/2016GL068009
An upper bound on the amplitude of mooring motion is estimated as follows. We decompose the observed velocities into mean and perturbation velocities: where subscripts i and c denote inner and center moorings, respectively. Then The first term is a secular signal while the second is a high-frequency signal representing the "flutter" of moorings. One can eliminate the former while retaining the latter by applying a high-pass filter to equation (A3). We used a cut-off frequency of 48 h, retaining all motions with periods smaller than this value. We next integrated this signal in time to obtain an upper bound on the zonal distance between the center and neighboring inner moorings: It is an upper bound since the mooring cable will constrain the magnitude of relative mooring motion to be smaller than that inferred from equation (A4). Finally, the variance, 2 ci x , is obtained by averaging x within a running window and estimating the deviation of x from this mean: where T=1 month. A similar calculation for v i and v c yields y and 2  Figure S4c illustrates the magnitude of relative mooring motion as a function of time.

A4. Quantifying Uncertainty
The stochastic model presented above was used to estimate uncertainty in statistics. To accomplish this, we employed a Monte Carlo approach. We first expressed perturbation distances in terms of time series, (t). Time series of ≡ ( x , y ) were created by randomly sampling Gaussian distributions with zero means and standard deviations, , given above. We smoothed resulting time series using a 24 h filter to prevent abrupt mooring motion. Next, the vector position of two moorings (cf. equation (2)) was defined as dl = dx + (t), where |dx| denotes the fixed distance between the moorings.
Having obtained stochastic realizations of mooring position, was calculated. For each time series, the standard deviation and skewness were estimated. Repeating this process 1000 times yielded distributions for each statistic as functions of depth and season. Shaded regions in Figure 2 illustrate the breadth of these distributions at the 90% level.