## Abstract

The local circulation of submesoscale fronts and filaments can be partly understood through a horizontal momentum balance of Coriolis, a horizontal pressure gradient, and vertical diffusivity in a turbulent boundary layer, known as the turbulent thermal wind balance (TTW or T^{2}W). T^{2}W often reproduces the instantaneous relative vorticity and divergence of submesoscale circulations in open-ocean and shelf settings. However, a diurnal cycle in submesoscale vorticity and divergence is characterized by a non-T^{2}W phasing: a maximum in divergence magnitude lags the maximum in vertical diffusivity (with vorticity lagging divergence). Here, an idealized model is used to solve the transient turbulent thermal wind (T^{3}W) equations and to investigate the diurnal evolution of front and filament circulation in a 2D plane. Relative to a steady-state circulation, transient evolution can cause both instantaneous strengthening and a weaker diurnal average of the secondary circulation. The primary mechanisms controlling the diurnal variability exist in a 1D Ekman layer that imprints onto the 2D circulation. In midlatitudes, acceleration during separate phases of the diurnal cycle (from night to day and from day to night) is dominated by distinct inertial oscillation and vertically diffusive mechanisms, respectively. However, the manifestation of these dual accelerations is sensitive to latitude, boundary layer depth, and the strength of the forcing. A simple 1D model predicts the diurnal phasing of submesoscale divergence and vorticity in realistic primitive equation simulations of the southwestern Pacific and coastal California.

## 1. Introduction

Submesoscale currents are relevant to a variety of oceanic processes (McWilliams 2016). These motions partially fill the spectral gap between the mesoscale and microscale to provide a forward route to energy dissipation (Barkan et al. 2015; Capet et al. 2008c; D’Asaro et al. 2011), influence density stratification in the surface layer (Boccaletti et al. 2007), affect the fate and transport of material in the ocean (Uchiyama et al. 2014; Gula et al. 2014; Romero et al. 2013), and exert control on ecosystem functioning (Mahadevan 2016; Lévy et al. 2012). These processes are controlled by the behavior of the recurrent spatial patterns (coherent structures) at the submesoscale: fronts, filaments, and vortices km in lateral scale that predominantly exist in the weakly stratified surface boundary layer (SBL) of the open ocean (Gula et al. 2014; Nagai et al. 2006; Capet et al. 2008b) and in the depth-filling, overlapping surface and bottom boundary layers of the continental shelf (Dauhajre et al. 2017, hereafter DMU17). The circulations associated with submesoscale structures are a hybrid of geostrophic and ageostrophic dynamics, quantitatively defined by Rossby and Froude numbers that are , where *V* represents a horizontal velocity; *f* is the Coriolis parameter; *L* and *H* are horizontal and vertical length scales, respectively; and is the stratification (where *ρ* is the density and *g* is the acceleration due to gravity). Physical investigation of submesoscale coherent structures has focused on regional phenomenology (Gula et al. 2014; Srinivasan et al. 2017; Bracco et al. 2016; Capet et al. 2008a; Kirincich 2016; DMU17), generation and destruction processes (Hoskins 1982; Boccaletti et al. 2007; Callies et al. 2016; McWilliams et al. 2015; Molemaker et al. 2015; Gula et al. 2014; McWilliams et al. 2009a; McWilliams and Molemaker 2011), seasonal sensitivities (Callies et al. 2015; Mensa et al. 2013), and controlling momentum balance relations (Nagai et al. 2006; Callies et al. 2016; McWilliams et al. 2015; McWilliams 2017).

Here, we extend the understanding of submesoscale coherent structures by exploring the mechanisms controlling a diurnal evolution in front and filament circulations. This diurnal variability has not been investigated in previous studies of submesoscale dynamics, yet it can impact the high- and low-frequency behavior of submesoscale circulations. These impacts on the circulation can in turn affect processes controlled by submesoscale flows described above. The problem is posed in an idealized setting and uses the turbulent thermal wind (TTW or T^{2}W) horizontal momentum balance relation (McWilliams et al. 2015; Wenegrat and McPhaden 2016b) as a starting point. T^{2}W is a diagnostic relation that assumes a balance among rotation, pressure gradient, and vertical mixing of momentum. Front and filament circulations can be calculated from this balance with knowledge of the transverse buoyancy gradient (here, defined as , where *b* is the buoyancy, the subscript denotes a partial derivative, and *y* denotes a cross-front axis) and parameterized vertical diffusivity profiles. For a filament, such a circulation is defined by a primarily geostrophic, cyclonically sheared longitudinal (along filament) flow and a transverse (cross filament) ageostrophic secondary circulation characterized by strong surface convergence and downwelling (Fig. 2; detailed in section 2a).

Relative to other frontal secondary circulation diagnostics (Garrett and Loder 1981; Nagai et al. 2006), T^{2}W is the most successful in describing the circulations of fronts or filaments with . A T^{2}W prediction for dense filament frontogenesis by boundary layer turbulence (McWilliams et al. 2015) compares well with a large-eddy simulation (LES) of the same process (Sullivan and McWilliams 2018). Similarly, T^{2}W successfully reproduces the local circulation of individual fronts or filaments in some regional simulations [e.g., Gulf Stream cold filaments investigated in Gula et al. (2014)] when is strong and there is a well-defined .

However, in analysis of a Regional Oceanic Modeling System (ROMS) simulation of the Southern California Bight, DMU17 show that there is a high-frequency, diurnal evolution of submesoscale circulations with three distinct analyses. First, it is shown that individual fronts and filaments can exhibit strong secondary circulations at times of weak ; these circulations are not accurately predicted by T^{2}W (Fig. 11 of DMU17). Further, the coincidence of strong secondary circulation and weak is characteristic of a domain-wide diurnal cycle in horizontal divergence *δ* and vertical relative vorticity *ζ* at submesoscale magnitudes . This domain-wide diurnal evolution is indicated by both a spatial root-mean-square (RMS) time series of secondary circulation strength indicators (e.g., , where *w* is the vertical velocity; Fig. 12 of DMU17) and Lagrangian sampled time series of *δ* and *ζ*, where particles are advected into fronts and filaments (Fig. 15 of DMU17).

Here, we present another metric evidencing a diurnal oscillation in secondary circulations in a realistic, full primitive equation (PE) simulation. Figure 1 shows the diurnal evolution of vertical diffusivity, vertical velocity, and relative vorticity local to a dense filament on the Santa Monica Bay shelf over a 24-h period. The dense filament is followed in time, and the fields local to its circulation are isolated with 3-km cross sections that cut across the density gradient defining the filament (further analysis details can be found in DMU17). Solar heating drives a diurnal cycle in (Fig. 1c). Converse to the TTW prediction of strong secondary circulation () at times of large , the strongest secondary circulation of the filament (Fig. 1d, dark red) occurs at a minimum in , with (Fig. 1d, dark cyan) lagging .

The above two paragraphs summarize four metrics evidencing a diurnal oscillation of submesoscale circulations in a midlatitude PE simulation. In all four analyses, the same qualitative phasing is apparent: a maximum in coincident with a minimum in and with lagging . These four separate indications of a diurnal, non-T^{2}W evolution in secondary circulations motivate investigation of the mechanism(s) controlling the diurnal evolution, the high-frequency and diurnally averaged effects on the circulation by transient oscillations, and the implications of this diurnal variability on submesoscale controlled processes (e.g., vertical fluxes in the boundary layer).

In this study, the basic physics driving the diurnal variability in front or filament circulation are explained with a modified T^{2}W system that includes acceleration, which we call the transient turbulent thermal wind (T^{3}W) system. The dynamical target is a mid–life cycle submesoscale density front or filament (after frontogenetic birth and before frontolytic destruction) in the surface boundary layer with and where the Ekman number fluctuates in magnitude over a diurnal period (here, is the surface boundary layer depth). Our goal is to understand the mechanism(s) that cause temporal evolution of front or filament circulation forced by a diurnal variation in solar heating or wind stress, the extent of the disagreement between T^{2}W and T^{3}W solutions, and the effect of the transient evolution on the low-frequency strength of the secondary circulation. We are motivated by this simple question: Why is there such strong when is weak in fronts and filaments? Further, is the diurnal evolution of the secondary circulation distinct or universal for different forcing regimes or geographical settings? Our analysis will expand on the success of the T^{2}W system in describing submesoscale front and filament circulations with high Ro (Gula et al. 2014; McWilliams et al. 2015; Wenegrat and McPhaden 2016b; DMU17), with this study correcting an apparent predictive limit of the T^{2}W balance during periods of high-frequency temporal variability in .

This paper introduces the T^{3}W equations as a 2D system that captures the relevant spatial structure of a front or filament circulation and accurately predicts the phasing and strength of secondary circulations under diurnal surface forcing observed in Fig. 1 (sections 2 and 3). However, the fundamental mechanisms controlling the diurnal evolution of the 2D circulation are shown to operate in independent, 1D vertical columns responding to temporal changes in (section 4). As such, this study adds and overlaps with the literature on the upper-ocean response to diurnal fluctuations in vertical diffusivity (McWilliams et al. 2009b; Price et al. 1986; Wenegrat and McPhaden 2016a), but it represents a particular, unexplored case of 1D Ekman layer dynamics translating to coherent 2D circulation structures. Our analysis will reveal and quantify the important role of both inertial and diffusive mechanisms in the 1D system and expose when, why, and with what magnitude differences between T^{2}W and T^{3}W solutions arise. In section 5, we will show that the phasing of *δ* and *ζ* in fronts and filaments observed in two ROMS simulations (southwestern Pacific and coastal California) exhibits regional sensitivities that are predicted by the 1D model. We conclude with a summary and discussion of results in section 6.

## 2. The transient turbulent thermal wind system

### a. T^{2}W balance

We pose the transient turbulent thermal wind equations starting with the diagnostic, turbulent thermal wind momentum balance (McWilliams et al. 2015; Wenegrat and McPhaden 2016b; McWilliams 2017). T^{2}W is contained within the hydrostatic primitive equations and combines geostrophic, Ekman, and hydrostatic dynamics to give a diagnostic relation for front and filament circulation:

where *x* and *y* are the longitudinal and transverse axes, respectively. The horizontal velocity vector is , *f* is the Coriolis parameter, is the vertical diffusivity associated with boundary layer turbulence, and is the pressure normalized by a reference density .

For the problems presented in this study, we assume a rigid lid and a finite depth. These constraints require an assumption regarding the barotropic component of the pressure field in (1) that we detail in section 2b. Hydrostatic vertical momentum balance and specified density and vertical diffusivity profiles can be used in (1) to give an estimate of front or filament circulation .^{1}

Given a and , (1) is capable of producing the defining characteristics of a submesoscale flow structure, illustrated in Fig. 2. In this example, the buoyancy structure (dashed lines) shows the dense filament structure: lighter fluid at the surface, with the transverse buoyancy gradient confined to the SBL () where strong vertical diffusivity (not shown) limits stratification (implied in the illustration is an assumption of longitudinal uniformity ). The associated T^{2}W circulation is composed of a primarily geostrophic, vertically sheared longitudinal jet that gives rise to cyclonic vorticity at the filament center and an ageostrophic secondary circulation (black streamlines) with surface convergence and downwelling at the filament center.^{2}

It is important to note that (1a) can be solved in individual vertical columns. In this manner, the horizontal structure of filament (or front) circulation (e.g., left panel of Fig. 2) results from adjacent 1D T^{2}W solutions that are dependent on individual vertical profiles and . Two examples of 1D T^{2}W profiles are shown in the right panels of Fig. 2, illustrating day (top; ) and night (bottom; ) T^{2}W predictions.

The profiles at large Ek (Fig. 2, top right) are typical of a midlatitude dense filament . The longitudinal geostrophic flow is eastward and decays with depth by thermal wind balance. The ageostrophic flow due to is a clockwise Ekman spiral: is positive (negative) at the surface and negative (positive) at the base of the boundary layer. The total longitudinal flow is reduced in magnitude from the geostrophic component by the ageostrophic component that is driven by the strong vertical diffusivity. At lower Ek (weak ), T^{2}W reverts to near-thermal wind balance with the ageostrophic flows and *υ* much weaker. Because of the symmetric structure of the filament (or a front), *δ* and *ζ* at the filament center can be thought of as proportional to the magnitudes of the 1D velocities [i.e., , where is a 1D TTW profile]. We will rely on this assumption in sections 4d and 5 to translate 1D ageostrophic dynamics back to the 2D circulation.

McWilliams (2017) derives a scaling estimate for the strength of the (overturning) secondary circulation based on the diagnostic T^{2}W system:

where is the transverse buoyancy gradient with associated geostrophic flow . Either large , large , or small produces strong secondary circulation overturning cells. Equations (1) and (2) are accurate in predicting the general structure and strength of front or filament circulation in a quasi-steady balance. Yet, the primary motivation of this study is based on observations in realistic submesoscale simulations that contradict (2) in the diurnal phasing of the circulation: strong at weak (Fig. 1).

### b. T^{3}W governing equations

Here, we introduce time memory to (1) to form the 2D transient turbulent thermal wind system. T^{3}W is posed in a transverse-depth plane under the assumption of longitudinal uniformity along a front or filament (, which excludes frontal instabilities other than symmetric). T^{3}W builds on (1) by adding time derivatives to the momentum equations and accounting for buoyancy evolution by advection and vertical eddy diffusion. These are the simplest a priori additions to T^{2}W that allow time memory. The hydrostatic assumption can be justified in solutions that are scale anisotropic (*H*/*L* ≪ 1, where *H* and *L* are vertical and horizontal scales, respectively). In the solution shown in section 3, surface and bottom boundary layer depths and vertical eddy diffusivities for momentum and buoyancy are parameterized with a *K*-profile parameterization (KPP; Large et al. 1994; McWilliams et al. 2009b).^{3}

Momentum advection is not included for simplicity, and the buoyancy advection allows some frontogenesis to occur, although less than if momentum advection were retained. In highly frontogenetic situations, momentum advection does play a significant role (McWilliams 2017), but this is not the situation being addressed in this study. Nonlinear dynamics due to momentum advection can potentially cause variability in secondary circulation over a diurnal period [e.g., by the nonlinear Ekman transport effect (Wenegrat and Thomas 2017) or by a favorably aligned, oscillating wind stress (McWilliams 2017)]; however, we do not explore such effects in this study.

The 2D T^{3}W system is given by

Definitions are as in (1), with *u*, *υ*, and *w* velocities in the *x*, *y*, and *z* directions, respectively; *f* assumed constant in *y*; ; and .

Vertical boundary conditions for momentum are

and buoyancy

where represent surface and bottom stresses, respectively. Term is a thermal expansion coefficient; is a surface heat flux (W m^{−2}); and is a specific heat capacity. The bottom stress is computed as a bottom-layer quadratic drag law:

where is the velocity vector at the bottom; is von Kármán’s constant; is the height above the bottom (midpoint of the bottom grid cell; see online supplemental information); and is a roughness length. The vertical velocity is obtained from vertical integration of the continuity equation with a rigid-lid assumption.

A diagnostic barotropic pressure equation is used to compute and is derived by taking the divergence of the vertically averaged horizontal momentum equations:

The overbarred variables in (7) denote the barotropic fields. That is, for any variable :

where the superscript *H* indicates a vertical average, and variables with a tilde represent a calculated baroclinic component; for example,

With the rigid-lid assumption, from continuity, and (7) becomes a diagnostic relation for :

The total pressure is then obtained by correcting the baroclinic pressure by (8).^{4} Equations (3) and (10) are solved numerically, and details of the spatial–temporal discretization are given in the online supplemental information (model code can be found at https://github.com/ddauhajre/T3W_DIM_MODEL).

## 3. Diurnal evolution of a 2D dense filament

Here, we show the success of (3) in capturing the diurnal variability in a 2D front or filament circulation. We will also show that the primary diurnal behavior is retained in a simplified, 1D version of (3) that is explored in section 4 and validated against ROMS data in section 5. The reduction in dynamical complexity is justified with two solutions. One solution (FIL1) contains as much dynamical complexity as allowed in (3). The vertical diffusivity contains spatial structure, is parameterized by KPP (along with boundary layer depths), and is forced by a diurnal heat flux (where *T* = 1 day and ) and a constant wind stress . Also, in FIL1, the buoyancy field evolves due to vertical diffusivity of buoyancy and advection. A second solution (FIL2) employs the following simplifications: vertical diffusivity is spatially uniform and prescribed , there are no surface stresses or fluxes , and buoyancy advection is turned off .

FIL1 and FIL2 are simulated at midlatitude to allow a qualitative comparison with the observation from Southern California in Fig. 1. The initial condition creation and model spinup for case FIL1 is described in appendix A. To allow a consistent comparison between the cases, FIL2 is initialized with the buoyancy field of case FIL1 after a 4-day spinup period, during which buoyancy advection has altered the idealized . Similarly, the magnitude of in case FIL2 is taken as a spatial average of the parameterized in case FIL1. We show FIL1 2 days after the spinup procedure, where the solution exhibits approximate periodicity (there is a low-frequency trend in by the secondary circulation advection). FIL2 has no such trend.

After the 4-day spinup, the dense filament is within a submesoscale parameter space with a local mean and for FIL1, where the mean is taken over a diurnal period in the 2 km surrounding the filament center and in the upper half of the water column.

### a. FIL1 diurnal evolution

Figure 3 shows the evolution of , , and secondary circulation streamlines for FIL1 over a 24-h period. Strong and overlapping (gray line) and (orange line) at the filament center are present at hour 0 (Fig. 3, first column). At hour 9, surface heating increases stratification that causes to decrease in magnitude. Yet, reaches its maxima at hour 9 (Fig. 3, bottom). Lateral buoyancy advection associated with stronger at hour 9 sharpens at the surface and widens at the bottom (Fig. 3, first row). At hour 18, and separate and shrink toward the surface and bottom, respectively, as the heating continues. Term is strongest at hour 18, with some spatial deformation resulting from the evolution of the buoyancy . The effect of this deformation also imprints onto the secondary circulation as seen in and the streamlines of the secondary circulation. Domain-wide is much weaker at hour 18 relative to hour 9. Cooling begins after hour 18, and the fields evolve back to a similar state as in hour 0. The qualitative phasing of in the FIL1 solution agrees with the observations in DMU17 and Fig. 1.

### b. Diurnal behavior retained in 1D

Figure 4 shows a 2-day RMS time series of and for FIL1 (solid) and FIL2 (dashed). The RMS for any field *q* is denoted by and denotes the RMS normalized by its standard deviation in time. Terms are periodic at this stage in their respective simulations. In Fig. 4f, we show the T^{3}W (salmon) and T^{2}W (maroon) to compare instantaneous differences between the prognostic and diagnostic predictions of the secondary circulation. The diagnostic T^{2}W solution at each time step is calculated with the , and from the T^{3}W solution, where is used in (10) to compute .

The T^{2}W prediction (Fig. 4f, maroon curves) of the secondary circulation in both cases disagrees with the respective transient solution in its phasing. The diagnostic prediction primarily follows , as dictated by (2), even in FIL1, where is not constant in time. FIL2 has no modulation by advection (as well as constant , and spatially uniform ). Yet, the phasing of the RMS fields for FIL1 and FIL2 are extremely similar. In both cases, (Fig. 4c) peaks with a minimum in (Fig. 4a), and (Fig. 4b) lags . The secondary circulation strength, indicated by (Fig. 4e), is very similar in both cases, with a maximum around hour 8 each day coincident with a minimum in . Because is not changing in FIL2, the agreement between FIL1 and FIL2 indicates that ageostrophic accelerations are controlling diurnal phasing. Again, note the agreement in phasing between the simple 2D and 1D systems in Fig. 4 and the realization from a realistic PE simulation in Fig. 1. The transient oscillations in the secondary circulation are instantaneously stronger than the T^{2}W flow at times of weak mixing; however, the diurnal average (not shown) of the secondary circulation is weaker in a T^{3}W solution relative to T^{2}W. This averaged weakening is analogous to diurnal Ekman layer rectification (McWilliams et al. 2009b; Wenegrat and McPhaden 2016a), where the “effective” vertical eddy diffusivity in the diurnally averaged transient solution is reduced relative to a steady state.

## 4. 1D controlling mechanisms

### a. 1D T^{3}W governing equations

To isolate the fundamental mechanism controlling diurnal variability of front or filament circulation, we simplify (3) and isolate the ageostrophic flow. We decouple the momentum and buoyancy equations in (3) by assuming , subtract out the geostrophic component, and assume (and constant in time):

where now, .

The boundary condition enforces geostrophic shear at the surface and bottom of the surface boundary layer for *u* and a free-slip condition at the surface and bottom for *υ*:

where the geostrophic longitudinal shear . The problem can be analogously posed with a no-slip condition at the bottom (not shown), which produces the same general solution behavior. We note that this system has been studied in the atmospheric context (Zhang and Tan 2002; Van de Wiel et al. 2010) and more recently explored for the oceanic boundary layer (Wenegrat and McPhaden 2016a). Although, here, we use a free-slip boundary condition and geostrophic relation at the boundary, as opposed to a surface wind stress. We will also use a flow decomposition and nondimensional framework to explore the solution, as opposed to the analytical framework of Wenegrat and McPhaden (2016a). These 1D Ekman layer dynamics are shown to behave by dual inertial and diffusive processes. Ultimately, we relate the 1D dynamics back to 2D secondary circulations (sections 4d and 5).

We assume that varies periodically with a bulk change and define a mean vertical mixing . We also assume that the barotropic component of (11) is unforced and choose solutions with .

The full solution is decomposed into steady and periodic components:

where represents the steady solution, and represents the periodic fluctuation about the steady state. Velocity component is assumed to have a time and vertical average of zero. The steady solution is given by

with boundary conditions

with as free-slip boundary conditions at the surface and bottom.

#### Nondimensionalization

To understand the parameter regimes controlling the evolution of the circulation, we nondimensionalize (14) and (16). We choose scalings , , , , , , and . Here, *T* represents a diurnal period (1 day) and a geostrophic shear. We introduce a nondimensional vertical diffusivity , where and is a nondimensional parameter that expresses the bulk change in vertical diffusivity relative to the time mean.

Defining and choosing , (14) can be nondimensionalized:^{5}

with boundary conditions and at . Equation (17) represents a background Ekman layer system that sets and passes on vertical structure to .

The periodic problem (16) becomes

with analogous boundary conditions to (16). Here, is a ratio of the inertial frequency to the diurnal period , and represents a ratio of the range of mixing time scales relative to the diurnal period. Parameters and express the control by inertial or diffusive dynamics on accelerations.

The periodic solution is further decomposed:

where represents the periodic T^{2}W component, and represents the difference between the periodic T^{3}W and T^{2}W solutions. Velocity component evolves as

with the same boundary conditions imposed on .

The evolution of is given by

From this decomposition, we see the same inertial and diffusive parameters on the lhs of (20) [analogous to (18)]. Flow component changes due to the time tendency of the mixing [the second lhs term in (21)] with a baseline vertical structure imposed by the steady solution [rhs term in (21)].

For completeness, and for reference to metrics used in section 4b, evolution of the full-flow is given by

The nondimensional parameters and *k* will control the solution behavior, with and exerting the most control on the solution behavior in physically relevant parameter spaces (defined below). Equations (17), (18), (21), and (20) are solved numerically (detailed in the online supplemental information; model code: https://github.com/ddauhajre/T3W_ND_MODEL). We solve for the steady-flow first, then time step for the periodic transient component . At each time step, the periodic T^{2}W component and periodic difference component are computed from the periodic transient component . We initialize . All of the 1D solutions shown in the following subsections (Table 1) are run for 10 diurnal periods and exhibit periodicity.

### b. Base case

The nondimensional system [(17)–(21)] takes four realistic parameters as inputs: , and . We present a base case within a physically relevant parameter space for a submesoscale front or filament undergoing diurnal forcing , which gives the nondimensional parameters (Table 1). This solution illustrates relevant 1D solution behavior and highlights the dual inertial and diffusive mechanisms driving temporal variability.

Figure 5 shows the evolution for the base case solution. The full flow shows oscillations in *u* and *υ* characterized by a deepening and strengthening of the near-surface and bottom flows. The vertical structure of these flows is controlled by through the rhs term in (19). In the base case, the primary oscillation in *υ* leads *u*, with the strongest *υ* occurring near the minima in the nondimensional mixing at and the strongest *u* occurring at . The periodic transient components differ greatly from the periodic TTW component —an indication of the failure of the TTW system to capture the diurnal motion, especially at midday , when mixing is weakest. This difference is demonstrated by the large magnitude of . The structure of *u* and *υ* is in general agreement with the solution of Wenegrat and McPhaden (2016a). In section 4d, the phasing of *u* and *υ* in this solution will be translated back to *ζ* and *δ* (in a front or filament) to illustrate the predictive capability of the 1D model.

Inertial or diffusive terms can be responsible for . One way to quantify the relative importance of these terms is to look at their relative magnitudes in the momentum balance (as in Wenegrat and McPhaden 2016a). Here, we take a different approach and diagnose the dominance of either with a difference of normalized covariances, based on terms in (22) and defined for *x* and *y* components as follows:

where ; the overline in the denominators denotes a time mean; and denotes an RMS in *z*. and imply that inertial motions cause the acceleration terms , and and imply that the vertically diffusive terms cause .^{6} In the analysis that follows, we show the sum of these terms , which leads to the same interpretations as analysis of the individual time series.

We also quantify a normalized magnitude of the difference of from as

Both instantaneously large and a time mean over a diurnal period *T * indicate disagreement between and .

Figure 6 shows (Fig. 6a) and *D* (Fig. 6b) for the base case. As decreases to its minima at , , which implies that inertial terms cause the acceleration during the night-to-day decrease in mixing. Conversely, as increases back to its maxima , , indicating that diffusive terms cause the acceleration during the day-to-night increase in mixing. The curve (Fig. 6b) shows that there is a persistent nonzero deviation of the T^{2}W and T^{3}W solutions, with the largest difference occurring after the minima in mixing after the inertially caused acceleration of the flows.

The analysis of the base case shows that inertial- and diffusive-driven accelerations control the phasing of *u* and *υ* on either side of the diurnal cycle. The T^{2}W solution is unable to capture these accelerations and incorrectly predicts the diurnal evolution for this parameter regime.

### c. Parameter variations

The main controls on the nondimensional system are expressed in parameters and , which set the diffusive and inertial time scales, respectively. Parameter *k* will primarily control the amplitude of the response, and sets the vertical structure of the solution. We present two more solutions, P1 and P2, which increase and , respectively, relative to the base case (Table 1). P1 increases the range of diffusive time scales with , and P2 increases the inertial time scale with . These increases in and are relatively large and are used to clearly illustrate the isolated effect of each parameter.

In P1, the increase in results in a sharper oscillation in that is more confined to the surface and bottom layers (Fig. 7a). Here, there is still substantial disagreement between the periodic transient and periodic steady components, indicated by relatively large at midday (Figs. 7b,c). In P2, the increase in results in an oscillation that is phase shifted relative to the base case (positive peaks earlier; Fig. 7d). The periodic steady and periodic transient components agree more in their phasing, but the difference between the two still exhibits large magnitude (Figs. 7e,f).

The increase in in P1 causes the diffusive-driven acceleration to appear at the start of the diurnal period (Fig. 8a, green curve). There is also a more exaggerated inertial acceleration midday at (Fig. 8a, green curve). P2 exhibits predominantly inertial dynamics with a peak midday, followed by a short diffusive-driven acceleration at (Fig. 8a, orange curve). However, the large in P2 dictates mostly inertial dynamics. Both cases show large midday deviations between the transient and steady periodic components (Fig. 8b). In P1, the largest deviation occurs earlier relative to the base case (Fig. 8b, green curve), and in P2, the persistent inertial dynamics cause a larger and more sustained deviation throughout the diurnal period (Fig. 8b, orange curve).

### d. 1D translation to δ and ζ

The 1D nondimensional solutions explored above can ultimately be translated back to a prediction of surface divergence and vorticity phasing in a front or filament circulation with the underlying assumption of slowly changing geostrophic shear . Equation (22) gives the evolution of the ageostrophic flow in such a circulation, and a translation to a nondimensional divergence and vorticity can be made as follows: and . Here, the divergence is completely ageostrophic, and the factor of 2 is to emphasize that the divergence is simply the difference between two oppositely signed vertical columns on separate sides of a dense filament. The cyclonic vorticity of a filament is modified by the ageostrophic vorticity that modulates the geostrophic component of the vorticity , which we set to 1 for simplicity (i.e., the cyclonic vorticity is mostly geostrophic and modulated by the anticyclonic, ageostrophic vorticity evolution).

Figure 9 shows these translations for the base case, P1, and P2 solutions [based on the surface velocity ], which can be qualitatively compared to the metric from a ROMS simulation in Fig. 1. The base case displays peak *δ* at midday, with *ζ* lagging *δ*. In P2, *δ* peaks at midday, though less strongly than in the base case (note the different *y*-axis scalings), and *ζ* and *δ* are less phase lagged. In P3, *δ* is strongest earliest in the day, and *ζ* peaks midday. These results suggest a nuanced view of the diurnal evolution of *δ* and *ζ* in fronts and filaments, with phasing controlled by inertial and diffusive time scales set by latitude, boundary layer depth, and the strength of the diurnal forcing.

## 5. Validation of the 1D model with ROMS data

Here, we show the predictive capability of the 1D nondimensional system by using it to reproduce the phasing of *δ* and *ζ* observed in ROMS simulations of two different regions: southwestern Pacific (SWPAC; ) and coastal California near Pt. Conception (CCAL; ). Details of solution setups are in appendix B. A snapshot of surface relative vorticity for each solution (Fig. 10) shows the ubiquity of high Ro submesoscale structures in each solution. A diurnal cycle in , *δ*, and *ζ* for each solution is shown in Fig. 11 as a spatial RMS of surface and boundary layer for two diurnal periods. The RMS curves are generated from subregions in each of the simulations (Fig. 10, small black boxes) during times where multiple submesoscale fronts and filaments exist forced by primarily solar heating (with weak winds). There is a notable regional sensitivity in the diurnal phasing of *δ* and *ζ*. At lower latitudes (SWPAC; latitude ≈ −15°), *δ* peaks later in the day, and *ζ* and *δ* are less phase lagged (≈1-h lag). At midlatitudes (CCAL; latitude ≈ 34°), *δ* peaks midmorning, and there is a ≈2–4-h phase lag between *δ* and *ζ*.

As noted in the previous section, the only inputs needed for the nondimensional 1D system are , and . We take an average of and (for the times plotted) from the subdomains used to generate the curves in Fig. 11 and set a based on the first diurnal period in the RMS time series in Figs. 11a and 11c. These physical inputs and resulting nondimensional parameters are given in Table 2 for the two regions. We run the 1D nondimensional model [with a sinusoidal forcing , as in section 4]^{7} for the two regions and plot the translated (defined in section 4d) in Fig. 12 for two diurnal periods.

The 1D model is able to reproduce the general behavior for each region. At lower latitudes (SWPAC; Fig. 12b), *δ* and *ζ* have a very small phase separation and peak later in the diurnal period (coincident with weaker mixing, as in the ROMS data). In this case, diffusive dynamics are more dominant (i.e., controlled by ) and and in the 1D system evolve more in phase with each other. The CCAL solution (Fig. 12c) exhibits a larger phase separation between and , with peaking earlier, when the mixing has just reached its minima. Relative to SWPAC, this regime has a stronger expression of inertial dynamics that drive a larger phase separation between and . These results indicate that the simple 1D model is able to capture the primary processes controlling the diurnal phasing of submesoscale currents in a fully 3D PE solution and highlights a latitudinal sensitivity to the diurnal phasing of *δ* and *ζ* in submesoscale circulations.

The simple 1D model can guide observational campaigns aimed at sampling fronts and filaments at times of maximal convergence and vorticity. Figure 13 provides an approximate roadmap (and summary of general phasing trends) for *δ*, *ζ* phasing, based on the CCAL 1D solution for a range of latitudes (Fig. 13a) and (Fig. 13b) to show the effect on the phasing of for these important parameters. At lower latitudes or shallower boundary layers, the phase lag between *δ* and *ζ* will be small, and they will peak after or close to the minima in vertical mixing. At higher latitudes or deeper mixed layers, the phase lag between *δ* and *ζ* increases, and *δ* peaks earlier than the minima in mixing.

## 6. Summary and discussion

The T^{2}W diagnostic is capable of predicting the low-frequency spatial structure of 2D front or filament circulations, given a strong, nonchanging, boundary layer vertical diffusivity and lateral density gradient. However, large diurnal changes in vertical diffusivity cause fluctuations in the strength of front or filament circulations that are not captured by the instantaneous T^{2}W relation. In this paper, the T^{2}W relation is expanded to include time memory to give the T^{3}W system, posed here for a 2D transverse-depth plane. This simple 2D system is able to reproduce the observed diurnal phasing in a realistic PE simulation (Figs. 1, 11). Inherent in the 2D T^{3}W system are 1D Ekman layer dynamics that control the phasing of the flow over a diurnal period by dual inertial- and diffusive-driven accelerations in response to changes in (independent of surface forcing). The oscillations in the local 3D circulations can be attributed to these 1D mechanisms imprinting onto the front or filament circulation that is quasi-stationary in its spatial structure, as shown with an application of the 1D model to 3D ROMS solutions (section 5).

The transient T^{3}W solutions disagree with a T^{2}W diagnostic, and the extent of the discrepancies is controlled primarily by the magnitude of the change in over a diurnal period. The T^{3}W secondary circulation can be instantaneously stronger than its steady analog; however, transients driven by time variability in can weaken the diurnally averaged secondary circulation. These results hold for a density front (as well as the filament analyzed in section 3). A 1D nondimensional framework, posed with an assumption of stationarity in front or filament secondary circulation structure, exposes the primary mechanisms at play over a diurnal period. In midlatitudes, an inertially caused acceleration drives an increase in secondary circulation strength that leads to a large difference between T^{2}W and T^{3}W solutions as mixing decreases. Diffusive dynamics will primarily control the non-T^{2}W accelerations when mixing increases. The bulk change in vertical diffusivity, relative to its diurnal mean , controls the amplitude of the transient oscillation. The relative phasing of *δ* and *ζ* is primarily controlled by the inertial and range of diffusive time scales relative to the diurnal period ( and , respectively). At larger (smaller) , diffusive dynamics dominate, and *δ* and *ζ* are less phase lagged. The importance of these diffusive dynamics in 1D Ekman layer systems has not been discussed in previous studies (Wenegrat and McPhaden 2016a; Zhang and Tan 2002; Van de Wiel et al. 2010).

The T^{3}W equations lack momentum advection. In real fronts and filaments, the effects of momentum advection are not small, especially during frontogenesis (McWilliams 2017). Here, we have removed this source of nonlinearity for simplicity in interpretation of the behavior in question. Momentum advection evaluated with the 2D model solution (FIL1; where this term is neglected) show that it is smaller than the dominant retained terms in the linear balance (not shown), but not so small that it is an entirely accurate approximation. McWilliams (2017) shows that the main effect of nonlinear (ageostrophic) advection on a front or filament is a strengthening of the secondary circulation when the wind is oriented favorably (northwest or southeast for a filament and southeast for a front, for the transverse *y* axis defined here). It is possible that in the real ocean, the strengthening of the secondary circulation by momentum advection due to a favorable wind stress simultaneously occurs with T^{3}W oscillations due to changes in , though this remains an open question. However, the success of the 1D nondimensional system in predicting the diurnal phasing of submesoscale circulations in realistic, full PE simulations shows that linear 1D Ekman layer dynamics capture the fundamental processes relevant to diurnal variability when there is a well-defined diurnal signal in .

With the 1D model [(17)–(21)], all that is needed to predict the phasing of *δ* and *ζ* in fronts or filaments (subject to clean diurnal forcing) is , and . A priori estimates of these four parameters inputted into the 1D system can give an approximation of times of maximal *δ* and *ζ*. Such knowledge can guide both observational campaigns focused on high-frequency sampling of ephemeral submesoscale currents or pollution cleanup efforts on daily time scales. Similarly, the 1D model can provide an indirect measure of boundary layer turbulence with high-frequency measurements of .

The diurnal fluctuations of the secondary circulation are historically not considered in discussion of lateral or vertical fluxes induced by submesoscale coherent structures. Both the instantaneous and time-averaged modifications of the circulation have implications for fluxes (horizontal and vertical) associated with submesoscale structures. The role of submesoscale currents in regulating transport pathways of tracers is an open question, and the results of this work imply a more nuanced view of material transport on daily time scales by submesoscale circulations.

## Acknowledgments

We thank Alexander F. Shchepetkin for guidance on the numerical techniques, Roy Barkan and Jeroen Molemaker for discussion of the problem, Jacob Wenegrat and Roy Barkan for insightful comments on an earlier version of this manuscript, Kaushik Srinivasan for the SWPAC ROMS solution, and Cigdem Akan and Lionel Renault for the CCAL ROMS solution. This research is supported by the U.S. Office of Naval Research (Grant N000141410626) and the National Science Foundation (Grant OCE-1355970).

### APPENDIX A

#### 2D Dense Filament Idealization: Initial Condition and Spinup

##### a. Initial condition

We create and and feed them into the T^{2}W balance to give an initial condition for *u* and *υ* (and by continuity *w*). Following McWilliams (2017), a general 2D surface boundary layer filament can be idealized as

with a boundary layer profile designed to give a buoyancy extrema in the center,

For case FIL1, we prescribe the following parameters:

and set a shallow-water domain with

For such a filament, the boundary layer depth deepens by at the center of the filament with . We set . For all simulations, we use a time step and horizontal spacing .

An initial guess of is fed into a T^{2}W iteration with the idealized . Vertical diffusivity is idealized as in McWilliams et al. (2015):

where is a baseline amplitude for the diffusivity; is a background diffusivity in the interior; limits the vertical column maximum of ; and is a velocity logarithmic singularity regularization constant as .

We choose the following parameters:

##### b. Spinup

We spin up the solution to allow the KPP vertical diffusivities and boundary layers to adjust to a well-behaved state, as well as to allow any transient adjustments in other fields to take place. Diurnal cycling of the surface forcing (for case FIL1, ) begins once the solution has reached an equilibrated state. The spinup for case FIL1 is as follows:

Solve for initial

*u*,*υ*with (1) given and with and .Prescribe constant for days 1 and 2.

On day 3, linearly decrease to and linearly decrease to .

Hold and constant through day 4.

Diurnally vary beginning on day 5.

### APPENDIX B

#### ROMS Setup: SWPAC and CCAL

The SWPAC and CCAL simulations are made with ROMS (Shchepetkin and McWilliams 2005), a hydrostatic, primitive equation model that uses KPP (Large et al. 1994; McWilliams et al. 2009b) for vertical mixing. Submesoscale circulations are resolved in both SWPAC and CCAL with a nesting procedure. Both solutions contain 1-h instantaneous averaged output and are hindcasts of winter months, with SWPAC beginning on 19 July 2007 and CCAL on 20 November 2006.

The SWPAC solution nesting procedure is analogous to the solution in Srinivasan et al. (2017). An outer grid encompasses the entire Pacific Ocean with initial and boundary conditions given by SODA climatology. Three more nests are made at , , and . This innermost nest (Fig. 10a) contains 50 vertical levels and is used in section 5 to obtain the RMS curves in Fig. 11. Surface forcing is detailed in Srinivasan et al. (2017). We note that the SWPAC solution atmospheric forcing contains daily variability in the winds and an analytically prescribed diurnal heat flux. This combination of forcing gives a relatively clean hourly response to diurnal heating (Fig. 11a) that drives a well-defined diurnal response in the submesoscale circulation (Fig. 11b).

The CCAL solution (, 60 vertical levels; Fig. 10b) is a child nest of a U.S. West Coast solution set with an outer grid that is nested down to , and 0.1 km. The CCAL solution contains tides (prescribed in the outer parent grid) and atmospheric forcing (winds, surface heat, and freshwater fluxes) that is given by a 6-km Weather Research and Forecast (WRF) Model simulation. Compared to the SWPAC solution, the diurnal variability in CCAL is less idealized with hourly winds and heating given by WRF, though we restrict our analysis to times of clean forcing and response (Fig. 11c).

## REFERENCES

## Footnotes

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JPO-D-18-0143.s1.

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

Equation (1) is not a complete representation of a diagnostic balanced model for front or filament circulation. To satisfy balance, an omega equation is also needed that accounts for geostrophic tendency, advection, and the dynamical effect of interior stratification to correct the *w* calculated from (1) (McWilliams 2017).

^{2}

A warm (light) filament would have the opposite T^{2}W circulation (surface divergence and upwelling), although these features are short-lived and not as ubiquitous as their cold counterparts. Thus, we focus on a dense filament in this study.

^{3}

McWilliams et al. (2009b) and Durski et al. (2004) find that KPP is suitable for diurnal forcing and finite-depth application, respectively. Previous T^{2}W investigations utilize KPP in idealized (McWilliams et al. 2015; McWilliams 2017) and more realistic settings (Gula et al. 2014).

^{4}

If the T^{2}W relation is used to diagnose flows from a full primitive equation model [e.g., as in Gula et al. (2014) or DMU17], the full pressure field is known and can be used directly as an input into (1). However, if idealized T^{2}W solutions are being created from prescribed *b* and in finite depth with a rigid lid (as they are in this study), is taken from *b*, and is computed from a that is derived from a thermal wind relation with the prescribed *b*. In this sense, a T^{2}W initial condition can be created from an initially prescribed .

^{5}

The choice of is arbitrary, as it divides out in the nondimensionalization and only appears in the boundary condition for . Variation in does not change general solution behavior.

^{6}

and are a form of normalized covariance, but the individual terms in (23) subtracted from each other are not “standard” correlation coefficients (i.e., bounded by ). The metrics are bounded, but their bounds are dictated by the RMS. The numerator’s use of the product of separate vertical RMS terms produces similar results and interpretations to an RMS of the product of two terms , with only the magnitude varying between the two forms of the metric (not the temporal structure).

^{7}

The curves in Fig. 11 exhibit a steplike temporal shape. We use a sinusoidal forcing for the 1D model because is defined with a time average of zero. We note that while the temporal shape of can control some of the solution behavior, it does not greatly affect the general predictive capability of the 1D solutions.