## Abstract

The diurnal cycling of submesoscale circulations in vorticity, divergence, and strain is investigated using drifter data collected as part of the Lagrangian Submesoscale Experiment (LASER) experiment, which took place in the northern Gulf of Mexico during winter 2016, and ROMS simulations at different resolutions and degree of realism. The first observational evidence of a submesoscale diurnal cycle is presented. The cycling is detected in the LASER data during periods of weak winds, whereas the signal is obscured during strong wind events. Results from ROMS in the most realistic setup and in sensitivity runs with idealized wind patterns demonstrate that wind bursts disrupt the submesoscale diurnal cycle, independently of the time of day at which they happen. The observed and simulated submesoscale diurnal cycle supports the existence of a shift of approximately 1–3 h between the occurrence of divergence and vorticity maxima, broadly in agreement with theoretical predictions. The amplitude of the modeled signal, on the other hand, always underestimates the observed one, suggesting that even a horizontal resolution of 500 m is insufficient to capture the strength of the observed variability in submesoscale circulations. The paper also presents an evaluation of how well the diurnal cycle can be detected as function of the number of Lagrangian particles. If more than 2000 particle triplets are considered, the diurnal cycle is well captured, but for a number of triplets comparable to that of the LASER analysis, the reconstructed diurnal cycling displays high levels of noise both in the model and in the observations.

## 1. Introduction

Submesoscale circulations are consistently formed at the surface and bottom ocean boundary layers on spatial scales ranging from 100 m to a few kilometers. They have life cycles of few days and are characterized by Richardson (Ri) and Rossby (Ro) numbers of *O*(1) (Thomas and Ferrari 2008; McWilliams 2016). They play an important role in the energetic of the ocean system, by providing a bridge toward energy dissipation through loss of balance from the inverse energy cascading quasigeostrophic mesoscales (from tens to few hundreds of kilometers), and the forward energy cascading three-dimensional turbulent microscales (from centimeters to few tens of meters) (Capet et al. 2008; D’Asaro et al. 2011; Molemaker and McWilliams 2010; Barkan et al. 2015).

The formation and evolution of submesoscale fronts and filaments have been successfully described and understood in idealized configurations in terms of a turbulent thermal wind horizontal momentum balance relation (T^{2}W; Gula et al. 2014; McWilliams et al. 2015; Wenegrat and McPhaden 2016a) that assumes a balance between Coriolis force, pressure gradient, and vertical mixing. T^{2}W is one mechanism by which fronts are generated—the other two being mixed layer instability and strain-induced frontogenesis—and describes quite accurately the momentum balance after a front or filament has formed. Within this framework, a seasonal modulation of the submesoscale flow can be easily understood. From a theoretical standpoint, given a strong lateral buoyancy gradient, the dependence of submesoscale (Ro > 1) fronts and filaments on the vertical mixing coefficient *k*_{ν} is contained in the T^{2}W horizontal momentum balance relation. T^{2}W predicts that the strength of the ageostrophic circulations should be directly proportional to the vertical mixing coefficient. The greater the vertical mixing and the deeper the surface boundary layer, for a given lateral buoyancy gradient, the more potential energy can be released and the more energetic the submesoscale flow will be. In most oceanic regions the strength of submesoscale flows displays a seasonal cycle that follows that of the surface boundary layer, peaking in winter when the surface boundary layer is deepest (e.g., Callies et al. 2015; Mensa et al. 2013).

Alongside a seasonal cycling, Dauhajre et al. (2017) pointed out that submesoscale circulations are subjected to diurnal fluctuations while investigating their formation and relevance over the continental shelf of the Southern California Bight with a high-resolution numerical model. The diurnal cycling of dynamical quantities such as divergence, vertical velocity, and vorticity stems from diurnal variations in vertical mixing within the Ekman layer (Wenegrat and McPhaden 2016a), in turn forced by variations in solar heat flux and/or wind. Indeed, near the ocean surface vertical mixing increases in response to surface cooling and/or an increase in surface stress. However, Dauhajre et al. (2017) noted that off the Southern California Bight individual submesoscale fronts and filaments exhibit strong secondary circulations at times of weak *k*_{ν} and vice versa, according to a basinwide diurnal cycle that is contrary to the prediction of the T^{2}W model.

To account for such cycling, Dauhajre and McWilliams (2018) extended the T^{2}W horizontal momentum balance relation to include time memory by adding acceleration, and the time evolution of buoyancy by adding vertical eddy diffusion and advection. The resulting two-dimensional transverse-depth plane (the plane across the density gradient) system of equations has been named transient turbulent thermal wind (T^{3}W) balance. Using such a system they show that the mechanisms governing the diurnal cycle of submesoscale circulations operate on one-dimensional (1D) Ekman columns and that ageostrophic accelerations control the diurnal phasing. The vertical profiles of the ageostrophic component of the cross- and along-frontal flow exhibit a distinct vertical structure that responds to changes in vertical mixing whereby a decrease in mixing (which increases shear) increases the near-surface magnitude of the velocity. This variation in near-surface magnitude can be directly translated to changes in the cross-frontal shear of the flow (i.e., vorticity and divergence) due to the horizontal structure of submesoscale fronts or filaments. For example, in a dense filament, the divergence can be thought of as the difference between two equal and oppositely signed vertical profiles of the cross-frontal flow (which is entirely ageostrophic); a simultaneous increase in the near-surface magnitude of the cross-frontal flow on each side of the filament results in an increase in near-surface divergence. The accelerations are driven by inertial and diffusive dynamics operating in concert, in agreement with the analysis by Wenegrat and McPhaden (2016b).

Overall, the T^{3}W framework allows for predicting the diurnal evolution of near-surface divergence and vertical relative vorticity, their phasing and relative lag, that are controlled by the inertial and diffusive time scales of the flow, in turn a function of latitude, turbulent boundary layer depth, and strength of the diurnal forcing. It is worth noting that the T^{3}W mechanism has been found in model simulations that rely on the *K*-profile parameterization (KPP; Large et al. 1994) to represent boundary layer turbulence. Model outputs that resolved it explicitly (Sullivan and McWilliams 2018) agree with KPP simulation on frontogenesis. The submesoscale diurnal cycling is therefore expected to be a robust feature in ocean simulations if the resolution is sufficient, independently of the choice of boundary layer parameterization.

The local nature of submesoscale circulations and the associated high vertical velocities of the frontal and eddies structures (Koszalka et al. 2009) control lateral mixing up to several kilometers and vertical absolute dispersion within and across the mixed layer (Zhong and Bracco 2013; Poje et al. 2014; Choi et al. 2017; Zhong et al. 2017; Liu et al. 2018), with consequences extending to the functioning of the marine ecosystem (Lévy et al. 2012; Zhong et al. 2012; Mahadevan 2016). The submesoscale diurnal cycle may provide a modulation of mixing that is mostly unexplored to date. Work so far has focused on its theoretical understanding and on simulations with climatological forcing fields. Here we extend previous investigations by exploring the imprinting of the diurnal cycle on surface-confined Lagrangian tracers at the ocean surface in the northern Gulf of Mexico (GoM) using drifter data and a suite of model runs.

In the GoM the pervasiveness of submesoscale structures in offshore waters in winter found in ocean models run at kilometer-scale resolution (Luo et al. 2016; Barkan et al. 2017) has been confirmed by a large drifter experiment, the Lagrangian Submesoscale Experiment (LASER), that took place in January–February 2016 (D’Asaro et al. 2018) in the DeSoto Canyon area of the northern GoM (Fig. 1). In this work we use LASER data to explore if and when a submesoscale diurnal cycling is found in the ocean. Additionally, we use a set of simulations in different configurations and resolutions to investigate 1) under which conditions a model may capture the observed signal, 2) how heat and momentum fluxes that force the ocean surface modulate the submesoscale signal, and 3) how Lagrangian sampling affects the detection of the diurnal cycle.

In section 2, results from the T^{3}W theoretical model applied to the Gulf of Mexico are shown together with a preliminary analysis of the drifter observations as motivation for the work. Section 3 introduces model configurations and integrations, while modeling outcomes are discussed in section 4. A targeted data analysis from both models and observations is presented in section 5. Section 6 summarizes major findings.

## 2. Theoretical model and drifter observations

We introduce the analysis of the drifter observations by briefly discussing the results of the theoretical 1D, nondimensional T^{3}W model proposed by Dauhajre and McWilliams [2018, their Eqs. (17)–(21)]. The T^{3}W model accurately predicts the diurnal phasing of submesoscale circulation in realistic, primitive equation simulations under light wind conditions. This model prescribes uniform vertical mixing and solves for the time evolution of vertical profiles of ageostrophic velocities in a front or filament; these profiles are translated to vorticity and divergence as described in Dauhajre and McWilliams (2018). The 1D model takes four input parameters: maximum and minimum vertical mixing over a diurnal period, latitude, and averaged surface boundary layer depth over a diurnal period. The 1D T^{3}W model is implemented with the bulk parameters characterizing the DeSoto Canyon in the GoM (latitude 28°N) in winter, that is, maximum vertical viscosity *k*_{νMAX} = 0.052 m^{2} s^{−1}, minimum vertical viscosity *k*_{νMIN} = 0.006 m^{2} s^{−1}, and surface boundary layer depth (BLD) = 60 m. The *k*_{ν} diurnal profile used in the 1D model is idealized as sinusoidal. We have verified that the shape of the profile may at most induce a shift by one hour in the peak, but does not impact the overall behavior. The diurnal cycle at the ocean surface is evident in both lateral convergence and vorticity. In the convergence field the diurnal cycle peaks around midday, and this coincides with a minimum in vertical mixing, while the vorticity cycle lags by nearly 3 h (Fig. 1). Using a metric that quantifies the influence of inertial and/or diffusive regimes on diurnal phasing [the *C*_{x} + *C*_{y} metric defined in Eq. (23) of Dauhajre and McWilliams 2018], we also found that in the GoM the acceleration of the ageostrophic circulations are inertially controlled as vertical mixing decreases, followed by a shift to a diffusive regime as vertical mixing increases back to its maxima at night time (not shown). This behavior is generic at the latitude of the basin.

### a. Drifter observations and preliminary results

The LASER experiment (D’Asaro et al. 2018), performed during winter 2016 in the framework of the CARTHE (Consortium for Advanced Research on Transport of Hydrocarbon in the Environment) project, released a total of approximately 1000 drifters in the area of the DeSoto Canyon (Fig. 1), mostly during localized launches of hundreds of drifters in areas of approximately 10 km × 10 km. LASER provides the largest observational dataset to date to probe the existence of a submesoscale diurnal cycle. The CARTHE drifters are 85% biodegradable (Novelli et al. 2017) and use a Spot Trace GPS with an accuracy of about 7 m. The drifter properties in terms of slip velocity have been extensively studied in laboratory (Novelli et al. 2017) and open ocean experiments (Laxague et al. 2018), comparing them with a suite of other instruments. During LASER, all the launched drifters were drogued at 60 cm from the surface but, because of high winds and waves, 40% of them lost their drogues during the first 7 weeks of the experiment. Two different algorithms that detect drogue loss based on the properties of differential velocities of the drifters and GPS transmission rate have been applied to separate drogued and undrogued drifters (Haza et al. 2018; Barkan et al. 2019).

We introduce the LASER analysis focusing on the first 3 days of drifter trajectories after one of the localized launches, named Large Drifter Array or LDA. They provide motivation for the model investigation. More detailed results based on the complete dataset are then presented in section 5.

In this first analysis we used only drogued drifters identified through the method proposed by Haza et al. (2018). Tests on a subsample indicate an accuracy of 94%–100%, with detection occurring within a time range of 0.5–3 h from drogue loss for 85% of the drifters. In addition, the drifters were quality controlled and linearly interpolated to a regular 15-min interval, computing also average velocity (Haza et al. 2018). The LDA launch targeted a submesoscale cyclone, identified in high-resolution infrared images, and its surroundings, and released 326 drifters in a “radiator” pattern in the DeSoto Canyon on 7 February, during a period of relatively calm wind, ~8 m s^{−1} (this value corresponds to the average wind during LASER). The wind increased during the following days, reaching 15 m s^{−1} by the end of the second day after release. The analysis is performed using drifter triplets to reconstruct the dynamical properties of the flow, that is, vertical relative vorticity normalized by the Coriolis parameter [*ς*/*f* = 1/*f*(∂*υ*/∂*x* − ∂*u*/∂*y*)], lateral divergence, also normalized [*δ*/*f* = (∂*u*/∂*x* + ∂*υ*/∂*y*)/*f*], and horizontal strain, defined here as $S/f=[\u2061(\u2202u/\u2202x)2+\u2061(\u2202\upsilon /\u2202y)2+0.5\xd7\u2061(\u2202u/\u2202y+\u2202\upsilon /\u2202x)2]/f$. We do so by using the least squares method introduced by Molinari and Kirwan (1975) and Okubo and Ebbesmeyer (1976). Applications of this method can be found, for example, in Niiler et al. (1989), Paduan and Niiler (1990), and Righi and Strub (2001). Some triplets were removed based on the scale and shape of the triangle they formed, to avoid errors in the reconstruction, as explained in detail in section 4b.

The LDA launch provided a total of 132 original triplets of approximately 1000 m size, and the dynamical quantities have been calculated with a 15-min time step. The time series of RMS quantities during the first 3 days after release are shown in Fig. 2. While the analysis may suggest the existence of an oscillation with a period of about one day, that can be seen more clearly in the divergence and strain, the signal is not obvious and is noisy.

After the first three days, no obvious oscillation can be detected in the Lagrangian reconstructions considering original LASER triplets. By then, however, the sampling is very nonuniform, with many triplets tightly closed together inside the submesoscale cyclone and others stretched too far apart to be used for the least squares calculations.

To investigate the influence of data sampling and forcing issues on the estimate of the submesoscale diurnal cycle theorized by the T^{3}W model, we use a suite of numerical integrations of increasing realism and, for the most realistic realization, we consider a varying number of surface-confined Lagrangian tracers.

## 3. Model setup and simulation details

We use the Regional Ocean Modeling System (ROMS) configured over the GoM. ROMS solves the hydrostatic primitive equations with a third-order upstream biased advection scheme (Shchepetkin and McWilliams 2005) and the KPP parameterization for vertical mixing.

We consider two setups. The first configuration, characterized by climatological forcings, allows for isolating the submesoscale diurnal cycling and the role of horizontal resolution in its representation, and has been described in detail in Barkan et al. (2017) and Choi et al. (2017). We focus on integrations representative of February (winter) conditions, at horizontal resolution of 1500 m (or LR_{1}, low resolution) covering the entire Gulf of Mexico and nested into a run of the North Atlantic, 500 m (IR_{1} case for intermediate resolution) covering the northern portion of the Gulf and nested onto LR_{1}, and 150 m (HR_{1} case for high resolution) nested onto IR_{1} over the DeSoto Canyon region, where the 2010 Deepwater Horizon oil spill took place (Fig. 3). These runs are forced by Common Ocean Reference Experiment (CORE; Large and Yeager 2009) monthly heat fluxes, Hamburg Ocean Atmosphere Parameters and Fluxes from Satellite Data (HOAPS; Andersson et al. 2010) monthly freshwater atmospheric forcing fields, and daily varying, climatological winds build upon QuikSCAT data (Risien and Chelton 2008). An empirical diurnal cycle based on the shortwave flux is imposed. Furthermore, the Mississippi River system outflow is simulated based on monthly mean volume flux data (Dai and Trenberth 2002) in the LR_{1} run, and on daily volume fluxes reconstructed by USGS for the year 2010 in the IR_{1} and HR_{1} cases. No tidal forcing is included. Near-inertial oscillations are not captured due to the climatological wind field and the absence of tidal forcing.

The second configuration, characterized by increased realism in the forcing fields, allows for investigating the effects of wind and heat fluxes variability on the diurnal cycle, as well as the Lagrangian sampling. In this configuration the horizontal resolution is 500 m (IR_{2}), offline nested within a coarser resolution (1000 m) run, LR_{2}. The LR_{2} model domain covers the Gulf of Mexico north of 24°N and uses the Hybrid Coordinate Ocean Model–Navy Coupled Ocean Data Assimilation (HYCOM–NCODA) ocean prediction system as 6-hourly varying boundary conditions from January 2015 to December 2016, following a yearlong spinup. The 6-hourly heat and momentum fluxes and daily freshwater (evaporation and precipitation) fluxes from the ERA-Interim reanalysis (Dee et al. 2011) force LR_{2}, while the Mississippi River system discharge is imposed according to the USGS daily data. On this run, the IR_{2} case is nested over a domain similar to IR_{1}. Tidal forcing is neglected also in this case, while near-diurnal—given the Gulf of Mexico latitude—inertial oscillations can be excited due to resonance with the wind forcing. In winter, near-diurnal oscillations in this basin have been characterized as predominantly due to tidal forcing, while the wind-driven inertial contribution to the diurnal oscillations is diminished compared to summer due to a deeper mixed layer (Gough et al. 2016). This second configuration increases the realism compared to the previous by incorporating the observed variability in the diurnal cycling of winds and heat fluxes, and by partially including the inertial oscillation contribution. It also allows for investigating the role of interannual variability in the riverine inflow through the comparison of 2015, a year with a “standard” annual cycle in the Mississippi–Atchafalaya discharge (low discharge in winter, maximum discharge in late spring), with 2016, when the winter months were characterized by flooding conditions in conjunction with El Niño. Those solutions, however, do not contain a full realization of the wind-induced internal wave regime, due to the hydrostatic approximation, the lack of tidal forcing, limited time and space resolution of the momentum forcing, and the finite grid resolution [see also Beron-Vera and LaCasce (2016) for a quantification of the underestimation of the internal wave field in the GoM in the HYCOM–NCODA solutions used here as boundary conditions].

In the following, we analyze the diurnal cycling of Lagrangian dynamical quantities focusing on late January and February (winter, LASER period). Lagrangian particles are released near the ocean surface in the IR_{2} simulation, are advected by the lateral, two-dimensional velocity field, and are representative of surface drifters. The DeSoto Canyon region (within dashed blue lines) is seeded with 19 992 particles with an initial particle distance of 0.01° (DS, for DeSoto deployment). We further release 2240 particle triplets uniformly in the DS region with the initial distance between the vertices in each triangle being 1000 m. The Lagrangian tracers are confined near the ocean surface and are advected offline by interpolating on the particle positions the ROMS hourly averaged horizontal velocities using Larval Transport Lagrangian model (LTRANS) v.2b (North et al. 2011). We verified that at the resolutions considered in this work the hourly averaged fields introduce only a small error in the Lagrangian statistics.

## 4. The submesoscale diurnal cycle in the numerical simulations

The submesoscale diurnal cycle in the ROMS integrations is investigated considering its imprinting on vorticity, strain and divergence. We will show at the end of this section that in the model simulations these quantities can be computed with acceptable accuracy from Lagrangian tracer positions if the number of tracers per unit area is sufficiently large.

### a. The diurnal cycle in the Eulerian fields

We first focus on the more idealized configuration, that is, the runs described in Barkan et al. (2017), where an empirical diurnal cycle of shortwave radiation is imposed and the submesoscale diurnal cycling can be readily isolated. These runs allow us to cleanly quantify the resolution dependence of the statistics of interest. Figure 4 shows the evolution in the IR_{1} run of shortwave radiation, BLD and vertical mixing coefficient averaged over the boundary layer, together with the root-mean-square (RMS) of the Eulerian *ζ*/*f*, *δ*/*f*, and *S*/*f* calculated at a depth of 4 m, averaged over the area corresponding to the HR_{1} domain. Plots display the evolution of all quantities for 2 weeks in February and then zoomed over a 2-day portion to highlight the hourly evolution.

Parameter *k*_{ν} is negatively correlated with the shortwave radiation, with a sharp decrease that starts around 0700 EST to reach its minimum at midday, corresponding to the shortwave radiation maximum. The diurnal cycle in heating and cooling in turn modulates the BLD, which increases over night and rapidly decreases when the sun heats the ocean surface. A further small modulation to the variability of both *k*_{ν} and BLD, noticeable around day 12, is provided by the wind stress. Surface divergence, and therefore vertical velocity, responds to the changes in vertical mixing by increasing whenever the slope of the *k*_{ν} curve over time is negative, and vice versa, as described in Dauhajre and McWilliams (2018). The ROMS Eulerian RMS fields broadly agree with the extrapolated surface divergence and vorticity signals obtained with the 1D model showed in section 2. The peak in *δ* is followed by that in *ζ* with a delay of about 2.5 h. The relative phasing of *k*_{ν} and *δ*, and of *δ* and *ζ*, depends on the inertial time scale and on the ranges of diffusive time scales; therefore, it varies with latitude and with the BLD.

The Eulerian diurnal cycling on vorticity and divergence is visualized in Fig. 5, where snapshots are shown at approximately the (daily) maximum and minimum the following day in the IR_{1} and HR_{1} cases. The strength of the diurnal cycle varies with resolution because the strength of the submesoscale ageostrophic circulations is proportional to the model resolution. While the general shape and size of the submesoscale eddies, vorticity filaments and stronger fronts is comparable at 500- or 150-m resolution, further finescale structures are found in the 150-m case in the interior of, and between, the cyclonic submesoscale eddies (Barkan et al. 2017, their Fig. 4).

Figure 6 shows the evolution over five days of the RMS of *ζ*/*f*, *δ*/*f*, and *S*/*f* in the LR_{1}, IR_{1} and HR_{1} integrations calculated over the HR_{1} domain in all cases. In the LR_{1} run the diurnal cycle is barely distinguishable and is most evident in the HR_{1} case. The time-mean RMS of vorticity grows from 0.6 in LR_{1} to 0.9 IR_{1} and 1.4 in HR_{1}, a nearly constant increase of about 50% for a factor of 3 in horizontal resolution improvement. The RMS of divergence grows faster, with a mean value varying from 0.2 to 0.6 and finally 1.3 at the highest resolution adopted, where a sharper peak emerges. With regard to this dependence and to the subsequent comparison with the observational analysis, it is worth noting that the absence of air–sea coupling in the simulations may be causing an overestimation of the strengths of both mesoscale (Minobe et al. 2008; Seo et al. 2016; Renault et al. 2016, 2017) and submesoscale currents (Renault et al. 2018).

In the LR_{2} (not shown) and IR_{2} simulations (Fig. 7), the resolution dependence is comparable to that found in the first set of runs investigated. The RMS of vorticity grows from 0.60 ± 0.15 in the LR_{2} case, for which the horizontal resolution is 1000 m, to 0.73 ± 0.15 in the IR_{2} case for 2016 and 0.82 ± 0.16 in 2015 (the RMS and standard error are calculated based on the February month-long integrations). The variability of the cycling, on the other hand, is much greater and results from the variability in the momentum and heat fluxes used to force these runs.

In Fig. 7 both diurnal and interannual changes can be identified. At interannual scales there is an overall slightly weaker submesoscale field in February 2016 compared to the previous year that is reflected in the RMS values. The interannual differences are due to the massive winter flood recorded in December 2015–March 2016. The excess freshwater from the Mississippi River system caused a decrease in the depth of the upper turbulent boundary layer compared to the previous winter, suppressing in part submesoscale activity (Luo et al. 2016; Barkan et al. 2017).

The bottom panel in Fig. 7 illustrates the links between wind and surface radiation fields and diurnal cycling of submesoscales, suggesting that wind pulses disrupt the diurnal cycle the most. The cycling emerges more clearly during periods of calm wind conditions, as for example between 18 and 23 February 2016, and is disrupted following wind bursts. The role of the wind field is further quantified by the correlations between the daily maximum in vorticity *ς*/*f*_{Max} and the corresponding wind stress value, and between the amplitude of the diurnal cycle in vorticity (*ς*/*f*_{Max} − *ς*/*f*_{Min}) and the mean wind stress value over that day. The correlation coefficients over the 30 days for the curves showed in Fig. 7 are *r* = −0.60 and *r* = −0.61, respectively. If analogous correlations are calculated against the shortwave radiation, coefficients drop to *r* = −0.30 and *r* = −0.26. The difference in correlations further indicates that sudden wind changes disrupt the diurnal cycle more effectively than irregularities in heat fluxes.

To further examine the relationship between submesoscale diurnal cycle, shortwave radiation and wind bursts, we performed three sensitivity runs. In the first case we substituted the observed radiation fluxes with their climatological counterpart without changing wind forcing; the resulting cycling was indistinguishable from that presented in Fig. 7 (not shown). In the other two simulations we replaced the observed winds with a constant low value imposed from 29 January 2016, plus a Gaussian shaped peak at 0900 or 2100 local time 3 February (Fig. 8). When the wind is weak and steady, a stable daily cycle can be identified in all quantities. This is despite the BLD becoming shallower than in the IR_{2} run due to the persistent low winds. Once the wind starts to increase, the BLD deepens independently of the time of the day, and the RMS of *ζ*/*f*, *δ*/*f*, and *S*/*f* decreases. After the wind peaks, the BLD deepening stops, and the RMS of the quantities begins to recover. During the wind pulse, the RMS of *ζ*/*f*, *δ*/*f*, and *S*/*f* follows the wind change and reaches a minimum at about the time of the wind peak. Such behavior precludes the detection of a submesoscale diurnal cycle. After the wind burst, the submesoscale diurnal cycle reemerges. A video displaying the vorticity and lateral divergence evolution through the wind event is available in the online supplemental material to this paper.

The behavior is broadly consistent with the prediction of the T^{3}W model. The increase in mixing shifts the velocities back to a state of reduced vertical shear, which translates to weaker vorticity and divergence. However the 1D nondimensional model cannot formally simulate a wind burst because it assumes a constant diurnal average of vertical mixing (Dauhajre and McWilliams 2018). The evolution of the vorticity and divergence fields, shown in videos (supplemental material), also reveals an overall weakening of the strength of submesoscale eddies and fronts, suggesting that the 1D Ekman layer dynamics does not provide a complete description during wind bursts.

### b. The diurnal cycle in the modeled Lagrangian fields

In Fig. 9 the Eulerian values of vorticity, lateral divergence and strain are interpolated at the locations of the 19 992 evenly distributed particles that are released at 0600 local time in the DS region in the IR_{2} simulation. Within few hours, the diurnal cycle emerges. The irregularity found in the Eulerian cycling is duplicated in the Lagrangian RMS. The Lagrangian values are greater than the corresponding Eulerian ones and the amplitude of the diurnal cycle increases by about 30%, and more so in vorticity (see Fig. 9 versus Fig. 7). The Lagrangian amplification is due to the fact that surface-confined tracers preferentially sample regions with elevated negative divergence and positive vorticity, as shown already through the analysis of probability density functions (PDFs) and structure functions (Choi et al. 2017; Pearson et al. 2019). In the approximation of a periodic or close three-dimensional domain, Eulerian fields such as *ζ*/*f* or *δ*/*f* in each density layer have zero mean in the absence of large mesoscale circulations of a preferred sign. This is verified with very good approximation in the DS domain in February, as shown in Fig. 10. Given that flows with energetic submesoscale currents are characterized by skewed distributions favoring strong convergence/positive vorticity values, the end result is a larger number of grid points with small but positive divergence and negative vorticity values (Choi et al. 2017). Lagrangian particles converge preferentially in submesoscale currents with strong negative divergence and positive vorticity, and as many as ~70%–90% of particles are found in convergent and cyclonic structures within 2 days from their release.

To compare the model results with observations, we evaluate how much deterioration of the cycle signal occurs whenever lateral divergence and relative vorticity are reconstructed using only Lagrangian data, that is, are calculated from the particle locations, as done in the drifter case. To this end we consider the 2240 particle triplets uniformly seeded in the DS region. We then calculate the dynamical quantities of interest using the least squares method introduced by Molinari and Kirwan (1975) and Okubo and Ebbesmeyer (1976), as done for the LASER drifters. The precision at which this method allows for reconstructing dynamical quantities depends on the error in the available data and even more on the triplet configuration (Berta et al. 2016; Ohlmann et al. 2017). In the LASER case the position error associated with the GPS precision is on the order of 7 m, while in the modeled case, a small—but not zero—error is introduced because the Eulerian data are saved as hourly averages (Keating et al. 2011) and is further amplified by the LTRANS interpolation. The biggest errors in the reconstruction come from the triplet shape, that is, whenever the aspect ratio (AR) or the scale of the triplet is small (Ohlmann et al. 2017). The AR here is defined as the ratio between the height and the base for longest side of the triangle following Barkan et al. (2019). Only triplets that meet the condition AR ≥ 0.1 and for which the longest side ranges between 1000 and 3000 m are retained in the divergence, vorticity, and strain calculations, as done also in Ohlmann et al. (2017). The minimum length for the longest side assures that this is at least two model grid points; 1000 m is also the target distance between drifters belonging to the same triplet in the LASER deployment. A brief discussion on the sensitivity of our analysis to the choice of these thresholds can be found in the appendix.

Figure 10 compares the reconstruction of the dynamical quantities using the 2240 particle triplets with the “true” cycling (i.e., obtained from the Eulerian modeled fields interpolated to the triplet locations) in the mean and the RMS over five days. Two Eulerian curves are plotted: “Eulerian Avg.” is the arithmetic mean or RMS over the whole DS region, while “Eulerian Intp.” is the mean or RMS of the Eulerian values interpolated at the particle locations. We expect that the Lagrangian reconstruction can only capture Eulerian Intp. due to the preferential sampling by the particles discussed earlier.

The results in Fig. 10 show that the reconstruction is quite accurate for vorticity and strain, while tends to underestimate the mean of the Eulerian Intp. convergence, despite the large number of triplets considered. The underestimation is due to the fact that strong convergent filaments quickly and effectively stretch the triplets that sample them, and most of them have to be discarded from the reconstruction due to the AR constraint. Triplets in regions of elevated strain and vorticity, such as submesoscale eddies, on the other hand, are more likely to be retained in the calculation because are less effectively stretched. This is quantified in the right column of Fig. 10, where the PDFs of the dynamical variables are shown for the triplets retained in the reconstruction of the diurnal cycle over the first three days and in the case of the Eulerian interpolation at the position of all particles. The deterioration in the reconstruction of the skewness and tails of the Lagrangian distributions is greater for divergence than vorticity or strain (*y* axes are logarithmic in figure).

The dependence of the results on the number of initial triplets is further investigated by considering first the 1120 triplets released in the south half of the DS domain, then the 560 triplets initially in the southwest 1/4 portion, and finally the 140 triplets initially in the southwest 1/16 portion. As shown in Fig. 11, the ability to reconstruct the diurnal cycle from the Lagrangian positions with respect to the Eulerian interpolated values in the IR_{2} run does not vary significantly with the number of triplets, but if only 140 initial triplets are used the oscillation cannot be detected clearly in the RMS and the reconstruction becomes significantly noisier.

The diurnal cycle in the model with 140 triplets in Fig. 11 can be qualitatively compared with the observed one in Fig. 2, obtained with 132 original triplets launched in the LDA. The values of the dynamical properties are higher in the observations for all the dynamical properties, suggesting that despite the high resolution the model still underestimates the strength of the submesoscale circulations. The relative strength of the quantities, though, is similar, with divergence exhibiting the lowest values and vorticity the highest.

The model analysis presented so far helps explaining the noise level found in the observations and the fact that the diurnal signal cannot be easily detected using Lagrangian measurements (see Fig. 2). First, the number of triplets used in the reconstruction decreases in time because of the Lagrangian sampling and triplet stretching, with a consequent increase of the noise. Second, the strengthening of the wind during the LDA deployment after the first couple of days likely contributed to the obscuration of the cycle. We notice, though, that there are some differences between observation and model set up. The LDA results are obtained from original triplets launched within a submesoscale eddy, while the model results are based on a larger variety of triplet realizations. This implies that the LDA observations are less prone to the selective difference between divergence and vorticity shown in Fig. 10, but are more prone to depend on the occurrence of specific events that can dominate the dataset obscuring the cycle.

## 5. The diurnal cycle from chance triplets in model and observations

The results in section 4 highlight the challenges in identifying the diurnal cycle using time series of dynamical properties following a limited number (order of a hundred) of original triplets. The nature of Lagrangian sampling at the ocean surface favors convergent areas where triplets stretch, or submesoscale eddies where they accumulate close together. Triplets with a too large AR or too little separation cannot be used in the computation of dynamical properties, and the number of useful triplets decreases quickly in time.

To partially avoid this problem and maximize the number of triplets used in the statistics, in this Section we adopt a different strategy. Instead of following original triplets, we consider a period of time *T*, during which at each time step of 15 min for the LASER data or 1 h for the model we look for “chance” triplets of a given size, selecting only those that are almost equilateral. One potential concern when using chance triplets is that the positions and initial velocities may be correlated (LaCasce and Ohlmann 2003), but previous work indicated that the derived statistics are not systematically different from those from original triplets (Morel and Larceveque 1974; LaCasce and Bower 2000). The selected triplets are used to reconstruct the dynamical properties, and the results are binned according to the hour of the day, in order to investigate the average diurnal cycle during *T*.

The analysis in section 4 also shows that the occurrence of strong winds tends to obscure the diurnal cycle. Therefore, we perform a conditional averaging on the chance triplets, separating the triplets identified at each time in two ensembles, those above and below certain wind speed thresholds. Preliminary tests have been performed on the observational dataset, and the value of 7 m s^{−1} (just below the average wind speed of 8 m s^{−1} registered during LASER) has been chosen as one that clearly separates the two ensembles. We notice that the threshold value is likely to be dependent on the specific environmental parameters, and also on the specific realizations sampled by the drifters.

Chance triplets were first obtained for *T* corresponding to the month of February 2016, when the highest concentration of LASER data was available, chosen as any three particles found at a given time within a radius between 1000 ± 100 m and AR > 0.1. For each binned hour, we have 263 triplets on average (ranges from 63 to 2502) for the ROMS simulation (based on the case with 19 992 particles uniformly released) and 430 on average (ranges from 348 to 520) for the LASER drifters. Note that, in model, the output frequency is 1 h, against 15 min for the LASER data. Thus, even with a larger number of particles released in the IR_{2} domain, the number of available chance triplets within every hour is comparable to the observations.

Results for both model and observations are shown in Fig. 12. The RMS time series of the reconstructed dynamical properties over the month of February are first shown (Fig. 12a), together with the time series of the wind forcing (Fig. 12b), while the hourly reconstructed *ζ*/*f*, *δ*/*f*, and *S*/*f* are shown in Figs. 12c and 12d for periods of strong (>7 m s^{−1}) or weak (<7 m s^{−1}) wind, respectively. Comparing the amplitude of the mean daily cycle of the three dynamical quantities, it is evident that the diurnal cycle is obscured during episodes of strong wind.

In Figs. 12e and 12f, we show the hourly averaged daily cycle for the observations. The drifter results are noisier and have higher magnitude (more than a factor of 2) than the model ones. This is not surprising given that model results are limited by the resolution of 500 m and by hourly averaging, while the real ocean sampled by the drifters has a significant variability also at scales of less than 1000 m (Berta et al. 2016; Ohlmann et al. 2017).

The wind effect on the diurnal cycle is confirmed by the data. The diurnal signal is suppressed in presence of strong winds (Fig. 12e), while a cycle is present in the three dynamical quantities for weak winds (Fig. 12f). The daily maximum in the cycle occurs first in the divergence (around 13 h), followed by the strain and finally by vorticity (around 15 h), with a sequence and a time range compatible with the results of the ROMS model. The timing is also broadly in agreement with respect to the theoretical model (Fig. 1); differences (around one hour) are not surprising given the complexity of the ocean processes and the variability of the environmental parameters in the GoM. Noticeably, the peak in vorticity value reconstructed from the drifters is more than twice as large as in ROMS.

Figure 12 provides the first observational evidence of a submesoscale diurnal cycle, but the cycling is far less obvious than in the theoretical predictions. To further assess the robustness of the observational result, we modified the dataset, the wind threshold and the triplet parameters. These changes assure that the observational datasets are the least correlated but nonetheless robust in relation to the number of chance triplets available (many more triplets are indeed considered in this second case). We adopted the dataset and procedure introduced in Barkan et al. (2019). The LASER observations collected between 20 January and 20 February 2016 were separated in drogued and undrogued using more stringent criteria, and a local polynomial regression or LOWESS method (Cleveland and Devlin 1988) was used to reduce the GPS position noise with a 30-min smoothing window. This curve fitting approach ignores outliers, defined as data points outside the three standard deviations range and can reduce the differences in magnitude between the observed and modeled signal. Chance triplets were then constructed using drogued drifters that at each time where in the scale range (750–2250 m), with AR in the range (0.5–1.0). In this analysis, the scale range considered is significantly augmented, while the allowable AR values are more limited. Triplets were separated according to the wind conditions using two thresholds: greater than 11.2 m s^{−1} or less than 5.6 m s^{−1} for high and low wind conditions, respectively. Few instances with less than 100 triplets were removed from the calculation and weighting was applied to the mean quantities to account for the fact that some drifters occurred in many triplets.

The diurnal cycle in the dynamical quantities for this alternative LASER dataset and the model, using the new thresholds and triplet criteria, is shown in Fig. 13. The amplitude of all dynamical quantities is consistently lower for high wind conditions in both the observations and in the ROMS reconstruction, and the amplitude difference in vorticity between LASER data and model is reduced. In terms of diurnal cycling, the timing and relative delay among quantities, quantified using a harmonic fit to identify the peak, is smaller than in the Fig. 12, but continues to be broadly consistent with the theoretical and model predictions. It should be noted that while the RMS error in the LASER drifters is small, as the number of chance triplets is generally large (up to 20 000 and commonly few thousands), the variance between estimates for a local time is large on different days. This stems from the fact that the drifters sample relatively few features on each of those days (e.g., the LDA release targeted a submesoscale eddy). While there are differences between Figs. 12e,f and 13c,d with regards to amplitude and strength of the signal, due to the different datasets and triplet definition used in the analysis, range of AR considered and especially wind thresholds, they both support the existence of a submesoscale diurnal cycle.

Overall the results presented in this section point to the challenges associated with observing and modeling the ocean surface at submesoscale resolution. Even in the event of very large drifter releases such as LASER, the uncertainty resulting from sampling only a few features is large. They also suggest that in the ocean there might be various mechanisms influencing the wind response, other than the response to the wind burst in Fig. 8, for instance related to the response to sustained winds and interaction between winds and preexisting submesoscale features (Berta et al. 2020; Thomas et al. 2013)

## 6. Conclusions

Following previous theoretical work on the diurnal cycling of submesoscale circulations, we investigated whether such cycling is found in drifter data collected in winter 2016 in the northern Gulf of Mexico as part of the LASER experiment. With the goal of attributing the diurnal variability of the LASER data, we investigated the representation of the diurnal cycle of submesoscale circulations at the ocean surface in Eulerian and Lagrangian statistics in ROMS simulations of the region of interest. Comparing three integrations forced by climatological daily varying atmospheric forcing for the month of February at resolutions of 1500, 500, and 150 m we showed that Eulerian dynamical quantities, such as relative vorticity, lateral divergence, and strain undergo a diurnal cycle that increases for increasing resolution. Divergence is the most sensitive variable to model resolution of the three analyzed, and we found no evidence for convergence in the Eulerian statistics, in agreement with the analysis by Barkan et al. (2017).

We then increased the realism of the runs forcing ROMS at 500 m horizontal resolution with observed daily riverine discharge and 6-hourly winds and heat fluxes, and imposing boundary conditions derived by a data-assimilative hindcast. Surface-confined Lagrangian trajectories were advected in these runs. Given that Lagrangian tracers tend to concentrate in submesoscale features characterized by elevated vorticity, strain, and lateral convergence (Choi et al. 2017; Pearson et al. 2019), we found that the maximum values and the amplitude of the diurnal cycle are further amplified in Lagrangian measures. Additionally, we confirmed the role of freshwater forcing combined with the turbulent boundary layer depth in modulating the strength of the submesoscale field in this basin (Luo et al. 2016; Barkan et al. 2017) by showing that the record high riverine forcing experienced in winter 2016 weakened the strength of submesoscale circulations and their diurnal cycle compared to 2015.

We then evaluated how well the diurnal cycle can be detected as function of the number of Lagrangian particles released. In all cases, original triplets were used and followed for 5 days. If more than 2000 particle triplets are considered, the diurnal cycle is well captured, especially in relative vorticity and strain. Triplets trapped in eddies seldom get stretched below acceptable levels for the reconstruction, so that a high number of these triplets is retained; triplets sampling convergent fronts, on the other hand, have a greater chance of being eliminated from the reconstruction, with the end result that the mean of the convergence field is underestimated. Whenever 500 or fewer particles are released without targeting specific flow structures, they are more likely to sample convergent frontal structures than submesoscale eddies because of their overall area coverage, and the diurnal cycle is best recovered in divergence than in vorticity or strain; the divergence reconstruction, however, becomes noisy past the first couple of days. This modeling result could not be confirmed in the LASER data, given that the drifter releases targeted specific submesoscale structures.

For a number of modeled triplets comparable to that of the LASER analysis, the reconstructed diurnal cycling displayed comparable levels of noise in the model and observations, with better agreement between the two whenever more stringent criteria are applied in the selection of the drifter data. Noticeably, however, all variables in LASER analysis displayed a cycle stronger that found in the IR_{2} model, suggesting that a horizontal resolution of 500 m may not be sufficient to fully capture the strength of the submesoscale circulations and their variability.

Finally, we identified quasi-equilateral chance triplets at each time step, binned accordingly to their time during the day, and applied conditional statistics, separating the triplets occurring at weak and strong winds according to different empirical wind speed thresholds. For both model and drifter observations, and independently of the observational dataset considered, a root-mean-square diurnal cycle is observed for all dynamical properties in the case of weak winds, while for strong winds the signal is obscured. It is concluded that wind bursts disrupt the submesoscale diurnal cycle (Figs. 12, 13). This is independent of the time of day at which the wind burst happens, as highlighted through the idealized sensitivity runs. The cycle, when detected, suggests the existence of a shift of approximately 1.5–2 h between the occurrence of the divergence and vorticity maxima, in qualitative agreement with predictions of the T^{3}W theory. Quantitative disagreement should be expected given the lack of tidal forcing and the underestimation of near inertial oscillations in the simulations, and the lack of observations of mixed layer depth at the LASER drifters.

Future work will explore the implications of the diurnal cycling on planktonic systems and on air–sea fluxes, including those of oxygen and carbon.

## Acknowledgments

This research was made possible by a grant from The Gulf of Mexico Research Initiative through the CARTHE consortium. R.B. was supported also by the Office of Naval Research (N000141812697), National Science Foundation (OCE-1851397), and Israeli Science Foundation (1736/18). AB was supported also by the National Science Foundation (OCE- 1658174). Data are publicly available through the Gulf of Mexico Research Initiative Information and Data Cooperative (GRIIDC) (doi:10.7266/N7CN720V, doi:10.7266/N7PK0DK2, doi:10.7266/N7F18X4S, doi:10.7266/N75H7DQ5, doi:10.7266/N79885FW, doi:10.7266/n7-np9a-p268, doi:10.7266/N7W0940J).

### APPENDIX

#### Sensitivity of the Lagrangian Reconstruction

The error in the reconstruction of dynamical quantities inferred from drifter positions has been found to be sensitive to the shape and scale of the drifter clusters in previous work (e.g., Ohlmann et al. 2017; Barkan et al. 2019). Here we briefly explore such sensitivity in relation to the diurnal cycling. We examined the fidelity of the reconstructions using different aspect ratios (AR values, defined in the main text, section 4b) and with/without limitation on the size of the triangle of particles considered.

Figure A1 shows the reconstruction of lateral divergence using triplets selected using different AR but no specific constrain on the triangle shape. When the AR threshold is small, more extreme values are produced in the reconstruction and the variance (RMS) is overestimated; the mean, on the other hand, is mostly insensitive to the AR choice. Very high/low divergence values and unreal RMS peaks can be found even when the AR threshold is set to 0.5 and more than 96% of triplets are excluded from the calculation. Constraining the triangle scale is therefore necessary to achieve a realistic reconstruction. By excluding triangles that are too small, a large portion of the extreme values is removed (Fig. A2f), and the RMS is closer to the Eulerian Intp. (Fig. A2e). However, the reconstructed RMS is generally underestimating the Eulerian Intp. The reconstruction represents the average computed from all triangles, so the variance will be reduced when the quantities are calculated from very large triangles. This underestimation becomes less obvious when very large triangles are also excluded (Fig. A2h). Additionally, to limit this bias between reconstruction and Eulerian Intp. in the modeled data, it is useful to calculate the mean value of divergence, strain and vorticity for each triplet in the Eulerian field before calculating the RMS and PDF. This implies comparing one (mean, across the three particle locations) Eulerian value to the Lagrangian reconstruction, which indeed provides one value per triplet. In Figs. 10–13, both techniques (removing the large triangle and using the Eulerian mean value for each triplet) were applied.

## REFERENCES

*Tellus*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Geophys. Res. Lett.*

*Nat. Commun.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Amer. Stat. Assoc.*

*J. Hydrometeor.*

*Science*

*Proc. Natl. Acad. Sci. USA*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Quart. J. Roy. Meteor. Soc.*

*J. Geophys. Res. Oceans*

*J. Phys. Oceanogr.*

*J. Atmos. Oceanic Technol.*

*J. Phys. Oceanogr.*

*J. Geophys. Res.*

*J. Mar. Res.*

*J. Mar. Res.*

*Climate Dyn.*

*Rev. Geophys.*

*Geophys. Res. Lett.*

*Geophys. Res. Lett.*

*Elem. Sci. Anth.*

*Ocean Modell.*

*Annu. Rev. Mar. Sci.*

*Proc. Roy. Soc. London*

*J. Phys. Oceanogr.*

*Ocean Dyn.*

*Nature*

*J. Fluid Mech.*

*J. Phys. Oceanogr.*

*J. Atmos. Sci.*

*Deep-Sea Res. I*

*Monitoring and Modeling the Deepwater Horizon Oil Spill: A Record-Breaking Enterprise*, Y. Liu et al., Eds., Wiley, 217–226

*J. Atmos. Oceanic Technol.*

*Geophys. Res. Lett.*

*Deep-Sea Res. Oceanogr. Abstr.*

*J. Geophys. Res.*

*J. Phys. Oceanogr.*

*Proc. Natl. Acad. Sci. USA*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Sci. Rep.*

*J. Mar. Syst.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Ocean Modell.*

*J. Fluid Mech.*

*J. Phys. Oceanogr.*

*Ocean Modeling in an Eddying Regime*, M. W. Hecht and H. Hasumi, Eds., John Wiley and Sons, 17–38, https://doi.org/10.1029/177GM04

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Geophys. Res. Oceans*

*Limnol. Oceanogr. Fluids Environ.*

*Sci. Rep.*