Reduced-Order Quasilinear Model of Ocean Boundary-Layer Turbulence
Joseph Skitka J. B. Marston Baylor Fox-Kemper




Joseph Skitka and J. B. Marston
Brown University, Department of Physics 

Baylor Fox-Kemper
Brown University, Department of Earth, Environmental, and Planetary Science











The combined effectiveness of model reduction and the
quasilinear approximation for the reproduction of the low-order statistics of oceanic surface
boundary-layer turbulence is investigated. Idealized horizontally homogeneous problems 
of surface-forced thermal convection and Langmuir turbulence are
studied in detail.  Model reduction is achieved with a Galerkin projection of the governing equations onto an
subset of modes determined by proper orthogonal decomposition. 
For less than 0.2% of the modes retained, the reduced quasilinear model is able to
reproduce vertical profiles of horizontal mean fields as well as certain energetically important
second-order turbulent transport statistics and energies to within 30% error.  
For intermediate sizes of the basis truncation some statistics
approach those found in the fully nonlinear simulations.  Thus basis reduction 
can actually improve upon the accuracy of quasilinear dynamics.  A test
model with a small total number of modes demonstrates the non-monotonic convergence 
toward the correct statistics as the size of the basis is increased.

 Reduced-Order Quasilinear Model of Ocean Boundary-Layer Turbulence
    Andrea Dal Corso
    December 30, 2023
===================================================================













§ INTRODUCTION



Reliable weather and climate modeling requires the accurate simulation of two turbulent components of the Earth system that have thus far only seen incremental improvements: the oceanic and atmospheric boundary layers (BLs).  Given that long duration turbulence-resolving simulation is centuries beyond the reach of present-day computers <cit.>,
large-scale numerical simulations of the Earth's oceans and atmospheres employ the hydrostatic primitive equations, 
in which vertical accelerations of fluid parcels are presumed negligible, 
enabling significant computational savings.  Unfortunately, predictions of vertical transport and mixing are generally unreliable as these depend crucially on smaller-scale, non-hydrostatic 3D processes – particularly turbulence within the planetary BL <cit.> – that are  neglected or phenomenologically parameterized.  Errors in the simulation of low clouds, for example, are believed to be the single most important source of uncertainty in climate models and have stymied effort at improvement <cit.>.  Similarly, wind- and surface-wave-driven Langmuir turbulence in the ocean surface BL should be incorporated into climate models for improved accuracy and robustness <cit.>, yet a recent ocean BL simulation with horizontal grid spacing near 4 meters required millions of supercomputer hours for only weeks of simulated time to capture a region only 20km × 20km <cit.>.

Given that convective turbulence mediates vertical exchanges both within and between the atmospheric and oceanic BLs, it is perhaps not surprising that general circulation models are very sensitive to the chosen parameterizations of vertical transport <cit.>.  The present unsatisfactory state-of-affairs is summarized in a review of numerical weather prediction, which asserts “the uncertainties inherent to physical parameterization, either from incomplete process understanding or the dilemma of representing the impact of unresolved processes on the resolved scales, may require a fundamentally different approach.  Elements of parameterizations or entire schemes are likely to require components that appear statistical to the large scales.
Examples are stochastic sampling of parameter probability distribution functions, stochastically driven sub-cell models, or super-parameterizations through embedding entire convection-resolving simulations at sub-grid scale …
How radical this approach needs to be is currently not clear” <cit.>.  

The systematic development of Direct Statistical Simulation (DSS) approaches to the oceanic boundary layer is our central objective here.  As a first step in that direction we employ the quasi-linear (QL) approximation <cit.> for which a non-local second-order closure (described below) is exact.  <cit.> investigated the QL approximation in a related model of the dry convective atmospheric boundary layer, demonstrating that QL can capture BL deepening and, with less accuracy, the turbulent transport of kinetic energy.   Two ocean-surface boundary layer models will be used for case studies: Surface-forced thermal convection and Langmuir turbulence <cit.>.    Both problems are statistically homogeneous in the horizontal, while Langmuir and, if subjected to wind shear, convective turbulence additionally are anisotropic in the horizontal <cit.>.   We emphasize that, owing to horizontal isotropy, this configuration constitutes a particularly challenging application for DSS, yet one of fundamental geoscientific importance.

Eschewing the traditional route of accumulating statistics from numerical simulation time series, DSS is a set of approaches for directly obtaining the statistics of dynamical systems rather than by their accumulation from numerical simulation 
<cit.>.   One form of DSS focuses on the equations of motion (EOM) for the equal-time cumulants, which are truncated at low order.  Consider modeling the evolution of a state vector q(r⃗, t) specified by a system of master partial differential equations (PDEs). For simplicity, further consider the case where the nonlinearities are quadratic so that the system can be written as

    q_t =  L[ q]+ Q[ q, q],

where L is a linear vector differential operator and Q is the operator that includes the nonlinear (quadratic) interactions such as those arising in the material derivative.  Specializing to models that are translationally invariant in one direction, henceforth denoted the zonal (or streamwise) direction, a standard Reynolds decomposition of the state vector into its zonal mean and fluctuation (denoted with an overbar and prime, respectively) is made by setting q =  q+ q^'.  
The average here is chosen to satisfy the Reynolds rules of averaging so that q= q   and   q^'=0.  
Averaging Eqn. (<ref>), using the linearity of L and assuming the averaging operator commutes with L, yields

    q_t    =    L[ q]+ Q[ q,  q] =  L[ q]+ Q[ q,  q].

Of course, the difficulty now is to describe the average of the nonlinear term in the equation for the mean (the first cumulant). At this point local closure schemes are often adopted <cit.>.  For example, in large-eddy simulations (LES) of turbulent flows,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      the average of this nonlinear term is generally parameterized using Fickian diffusive schemes in which fluctuation correlations are taken to be proportional to the (local) gradient of mean variables <cit.>.  Not only are such closures inaccurate and not mathematically justifiable for inhomogeneous and anisotropic flows, but in LES the ad hoc proportionality coefficients must be trained for particular flow regimes.  Standard LES schemes are particularly ineffective near walls (when additional ad hoc, e.g. Van Driest, wall models are required) and for strongly anisotropic flows, including strongly stratified turbulence. 

A more sophisticated approach is to derive an equation for the evolution of the mean of the non-local fluctuation–fluctuation (or eddy–eddy) interactions <cit.>.  Unwarranted assumptions about homogeneity or isotropy are avoided, and non-local correlations of the form q'(r⃗_1) q'(r⃗_2) are tracked, where r⃗_1 and r⃗_2 are two position vectors within the domain of interest <cit.>. The equation governing the evolution of this quantity (the second cumulant) is obtained by multiplying the EOM  for the system (defined at the point r⃗_1) by q(r⃗_2), averaging and symmetrizing.  In addition to introducing non-local quadratic terms from the q  L[ q] terms, this procedure requires the evaluation of non-local cubic terms (the third cumulant) that arise from the q  Q[ q,  q] terms in Eqn. (<ref>).  The simplest and most computationally efficient approach is to truncate the hierarchy at second order by dropping the contribution of the third cumulant, q'(r⃗_1) q'(r⃗_2) q'(r⃗_3), to the EOM for the second.  Termed CE2 in <cit.> and S3T in <cit.>, this truncation of the cumulant expansion (CE) at second order is both a conservative and realizable approximation <cit.>.  Furthermore, CE2 is quasilinear in the sense that it includes interactions of mean quantities with eddies to give eddies and interactions of eddies with eddies to give mean flows, but neglects eddy-eddy scattering.  Equivalently, it is readily shown that by invoking a QL approximation of the instantaneous dynamics, 
the same closed equations for the first two cumulants are obtained <cit.>. 

While ineffective for homogeneous isotropic turbulence, the quasilinear approximation has been shown to qualitatively (and in some cases, quantitatively) reproduce mean flows and two-point correlation functions for model problems describing a wide range of anisotropic and heterogeneous flows.  Specifically, QL has been successfully utilized to describe the driving of mean flows in plasmas and on giant planets <cit.>, the sustenance of wall-bounded shear-flow turbulence <cit.>, the growth of a dry atmospheric convective boundary layer <cit.>, and even the development of the magnetorotational instability in accretion discs <cit.>.  The QL approximation is thought to be suitable to planetary flows because of the relative importance of interactions between small-scale
eddies and mean-field flow features <cit.>.  Recently, the suitability of this approximation for planetary boundary-layer
applications has been supported by demonstrations of the emergence of zonal jets on a randomly or stochastically forced β-plane <cit.> and the statistics of thermal instabilities in the atmospheric boundary layer <cit.> 
under QL dynamics or statistically equivalent frameworks.  Furthermore, the QL approximation has substantial capacity for improvement through
variations on the exact statistical closure of the cumulant expansion of QL systems (such as CE2.5 and CE3 <cit.>) and extensions such as the generalized quasilinear approximation (GQL) <cit.>.  

The quasilinear approximation has a number of mathematically attractive attributes.


  * QL can be shown to be exact in appropriate asymptotic limits (e.g. for strong mean flows or for an extreme separation of time scales) <cit.> and, in regimes away from the strict asymptotic limit, errors can be controlled and reduced by methodically retaining higher-order cumulants.

  * QL retains the quadratic conservation properties of the underlying master PDEs and is realizable, implying a non-negative probability distribution.

  * The low-order statistics found with QL vary more smoothly than do instantaneous flow-field realizations, suggesting the possibility of dimension reduction in numerical simulations <cit.>.

  * The attractor for the low-order statistical dynamical evolution may be much simpler (e.g. an equilibrium point or limit cycle) than that characterizing the instantaneous (generally turbulent) dynamics, and the temporal evolution of the low-order statistics is comparably stiff, implying that these simple attractors should be mored readily accessible. 

For all its virtues, the quasilinear approximation nevertheless suffers from two primary deficiencies.


  * Insufficient accuracy:  QL is deficient as the system is driven harder, reducing time-scale separation between the strict mean flow and the faster dynamics. This limitation was demonstrated in <cit.>, which compared the statistics derived from DNS for the problem of driving β-plane jets.  QL reproduced both the number and strength of jets for large time-scale separation but failed as the system was driven further away from equilibrium.  In practice, attempts to correct this deficiency by retaining higher-order cumulants are fraught with difficulty (e.g. the non-realizability of CE3).
More generally, higher-order cumulant expansions rapidly become analytically and computationally unwieldy.

  * Computational expense:  Although infrequently discussed in the research literature, naive implementations of QL/CE2 schemes actually can be as or even more expensive than full DNS because the second cumulant may have higher dimensionality than the underlying fields (depending upon the symmetries of the problem and the choice of averaging operation).  Extending the level of truncation by one order requires the calculation of correlations with an extra point in space, increasing the degrees of freedom of the system by D where D is the dimensionality of the system.  Consequently, in practice, attempts to increase the accuracy of CE2 by retaining higher-order cumulants are computationally prohibitive.


We show below that both of these deficiencies can be addressed simultaneously by  dimensional reduction of the boundary-layer models using Proper Orthogonal Decomposition (POD). POD, a method closely related to Empirical Orthogonal Functions (EOF) and Principal Component Analysis (PCA), is a common means to generate a reduced basis for applications in computational fluid dynamics in which the representation of a flow is decomposed into a set of modes based on the amount of variation in the flow that they describe in terms of energy <cit.>. Often, POD reduced models will retain only the most energetic of these modes. When the QL approximation is defined with respect to a spatial average over statistically homogeneous dimensions, it eliminates the direct flow of energy between scales in those dimensions; as a result, energy tends to populate fewer modes than with fully nonlinear (NL) dynamics. The use of energy optimization is a particularly suitable pairing for the QL approximation then because fewer modes must be retained to capture a similar portion of the variation in the flow than with NL dynamics.  

Another motivation for combining POD with QL when defined with a spatial average over a homogeneous dimension is that the statistics are invariant under a phase shift of individual Fourier modes in that dimension, implying that the phases of different Fourier modes are uncorrelated and localized structures cannot be captured in the flow solutions (see section <ref> of this paper, also <cit.>).  Because POD modes must be horizontal Fourier modes due to translational symmetry in the ensemble of flow realizations (see appendix <ref>), POD washes away localized coherent structures that could otherwise be represented with very few coherent modes.  (<cit.> notes a 𝒪(10^3) - 𝒪(10^4) reduction in the POD dimensionality of Rayleigh-Bénard convection when rigid walls break the symmetry of an otherwise periodic case.) The resulting bases have no encoded phase correlations, making the application of POD to homogeneous dimensions suited to the application of quasilinear approximation to spatial averages. 

Recently, <cit.> have demonstrated that accurate statistics for 2D rotating geostrophic flow on a sphere (a zonally-homogeneous, anisotropic, stationary problem) can be retained at massive dimensional reduction with a POD reduced model of CE2, an exact statistical closure of QL. Additionally, NL POD reduced models often suffer from a buildup of energy because energy pathways to the myriad small scales must be omitted, similar to an LES-simulation near the cutoff length <cit.>. In QL POD reduced models, small scales only exchange energy with with the mean flow along statistically homogeneous dimensions, thus there can only be an indirect buildup of energy at cutoff scales through the modification of the mean flow rather than a direct local buildup between scales. In light of this, it might be expected that QL POD reduced models would have better stability properties than NL POD reduced models. 

The particular problems used as test cases here, surface-driven
convection into stratification and Langmuir turbulence, are classic
problems in planetary boundary layers
<cit.>.  They are of immediate
importance to global climate modeling in both the ocean
<cit.> and atmosphere
<cit.>. From the oceanographic perspective, other
problems such as submesoscale mixed layer eddy restratification
<cit.> or mixing by symmetric instability
<cit.> are seemingly similar: inhabiting
the boundary layer and affecting the balance between mixing and
restratification. However, from the generalized symmetry perspective
submesoscale turbulence differs fundamentally, as it only exists in
the presence of horizontal density gradients, thus these problems have
horizontal homogeneity in only one direction.  The general approach
here applies to boundary layer turbulence that oceanographers call
“vertical mixing” and atmospheric modelers call “pencil models”, both of which are horizontally homogeneous in both directions and
heterogeneous in the vertical.

In one important way, the approach taken here resembles that of
<cit.>: it relies on a structured sequence of approximations
of increasingly higher moments.  However, Mellor-Yamada SGS
models <cit.> tend to emphasize moments that are
co-located in vertical position.  Doing so reduces the dimensionality
of the problem, but it prevents non-local effects from playing a
leading role in the scheme dynamics.  Non-local effects and relative
motions are fundamental to the statistics of turbulence, as frequently
demonstrated using structure functions and two-point correlations
<cit.>.  However,
allowing a SGS scheme to depend on non-local statistics means 
that the dimensionality of the analysis increases.  For example, for
a horizontally-homogeneous, stationary problem, one-point statistics 
will depend only on vertical position, while two-point, equal-time 
statistics depend on two independent vertical positions and their horizontal separation vector. Thus, a non-local approach is 
potentially more powerful, but not obviously efficient.

This article first describes the theory underpinning the proposed
model framework in section <ref>, including that of POD
Reduced Modeling, the QL approximation, and how basic assumptions
appropriate for ocean surface boundary layers, specifically horizontal
homogeneity, can be used to make simplifications.  Section
<ref> covers the numerical implementation of the this
framework, including an introduction to a custom-built reduced-model solver.
Section <ref> defines the two examples of ocean surface
boundary-layer turbulence that are used to test this technique.  The
results of the QL reduced-model runs are presented in section
<ref> while discussion of the implications of these
results as well as promising areas of future research are reserved for
section <ref>.








§ THEORETICAL FRAMEWORK



In POD, the flow is decomposed using the eigenvectors of the
uncentered, equal-time covariance matrix, 𝖱.  Let the flow
variables be written 𝐪 = (q_1(t)   ,    …    q_N(t))^⊤, where the index runs
over discretized spatial coordinates as well as fields (e.g. velocity,
temperature), t is the undiscretized time, and N is the total
number of flow variables.  
Returning to the Reynolds decomposition introduced in section <ref>),
𝐪 = 𝐪 + 𝐪' where
𝐪 indicates the mean fields under a generic
averaging procedure and 𝐪' indicates the fluctuating
component of the fields, then

    𝖱_ij = q_i q_j ; 𝖱 = 1/M∑_1≤ s≤ Mq_sq_s^⊤

where the indices i and j indicate both field and position.  s
indicates the field snapshot, of which there are M; these are
sampled from different times and/or ensemble members.  The
eigendecomposition of 𝖱 is:

    𝖱Φ^α = λ^αΦ^α

A truncated basis retains the most energetic modes,
Φ^α, as indicated by the associated
eigenvalues, λ^α.  The naked ϕ^α will be reserved for the
amplitudes of a flow field in the POD representation, whereas
Φ^α can be thought of as the unit vector of a
basis for flow states,

    q_i = ∑_1≤α≤ Nϕ^αΦ_i^α,
  
for some coefficients ϕ^α, where N ≤ M is the number of
independent modes spanned by the M snapshots. 

Because the cases of SGS turbulence are intended to be small in
horizontal extent, ignoring horizontal inhomogeneities in the flow
forcing is a reasonable approximation that will be important for
allowing the POD reduced model to be performed efficiently on a large,
three-dimensional problem.  For horizontally homogeneous cases, which
will be used exclusively in this work, covariances may be constrained
to depend on only the separation vector between two flow variables, 
 and it follows that POD
modes are constrained to be horizontal Fourier modes (see appendix
<ref> for a derivation of this.)  This allows for the
POD eigendecomposition (equation <ref>) to be performed
individually for each horizontal Fourier mode. To apply this to
discrete flow fields as they are represented in a solver, the fields
are written in terms of their horizontal Fourier wavenumber vector
index, 𝐦_i = {m_i,l_i} and their vertical
position and field index, z_i.  Here m_i and l_i are the
horizontal, x and y directed wavenumber indices.  Individual POD
modes from equation <ref> now need only be labeled with
their horizontal wavenumber and an additional index, k_α,
corresponding to the degrees of freedom associated with the vertical
direction and the different fields: Φ^α
= Φ_𝐦_α,k_α and may be defined by a
discrete vertical profile:

    Φ^α = (Φ_𝐦_α,k_α,z_1  ,    …    ,    Φ_𝐦_α,k_α,z_N_z)^⊤


Let the evolution of field variables, q, be written in terms of
generic constant, 𝖢, linear, 𝖫, and quadratic,
𝖰, operators in the horizontal Fourier basis and vertical
position basis.

    ∂_t q_𝐦_i,z_i(t)    = 𝖰_𝐦_i,z_i;𝐦_j,z_j;𝐦_l,z_l q_𝐦_j,z_j(t) q_𝐦_l,z_l(t) 
          + 𝖫_𝐦_i,z_i;𝐦_j,z_j q_𝐦_j,z_j(t) + 𝖢_𝐦_i,z_i

where t is the undiscretized time and repeated indices imply a sum.
Hereafter, t will be dropped and all field variables will be assumed
to have an implicit time coordinate and for the remainder of this
section, z will run over vertical positions and fields.  Note that
many elements in the implied summations will be exactly zero, allowing
for the equations to be constrained in terms of the allowable
horizontal wavenumbers in subsequent equations. 

As an example of the generic operator decomposition in equation
<ref>, consider a basic temperature advection equation,
∂_t T + (𝐮·∇) T =
F_T + κ∇^2 T, where T is the temperature, 𝐮
is the 3-dimensional velocity field, F_T accounts for external
forcing, and κ is the thermal diffusivity of a fluid. Here,
𝖰 would be a tensor strictly defined by the discretized
advective term, 𝖫 would be defined by the diffusive term,
and 𝖢 would be defined by the external forcing.
In fluids, the continuum equations that govern
fully resolved flows do not have terms that are higher than second-order
in the field variables, so this generalized notation is an accurate
starting point for developing a SGS model.

The quasilinear approximation 
neglects the scattering of eddy fields, q' =
q_𝐦0, into other eddy fields,
thereby eliminating the scale-by-scale transfer of energy between horizontal length
scales; rather, all energy transfer to or from a given length-scale of
turbulence must be with the mean field, q =
q_𝐦=0.  This omits the scale-by-scale cascade of turbulent
energy and enstrophy.  Assuming the
evolution coefficients 𝖰, 𝖫, and 𝖢
are also horizontally homogeneous, the quadratic and linear operators
are constrained to just the interacting “triads” of the quadratic term interactions within vertical water columns:

    𝖰_𝐦_i,z_i;𝐦_j,z_j;𝐦_l,z_l   = δ̃_𝐦_i - 𝐦_j - 𝐦_l𝖰_𝐦_i,z_i;𝐦_j,z_j;𝐦_l,z_l
    𝖫_𝐦_i,z_i;𝐦_j,z_j   = δ̃_𝐦_i - 𝐦_j𝖫_𝐦_i,z_i;𝐦_j,z_j

where δ̃_𝐦 is a modified Kronecker delta:

    δ̃_𝐦 = δ_MOD[𝐦_i],0 = δ_ m   mod   N_x,0δ_ l   mod   N_y,0

Then the subset of the coefficients not constrained to be zero within
this context may be written 

    𝖰_𝐦_i,𝐦_j,z_i,z_j,z_l   = 𝖰_𝐦_i,z_i;𝐦_j,z_j;𝐦_i-𝐦_j,z_l
    𝖫_𝐦_i,z_i,z_j    = 𝖫_𝐦_i,z_i;𝐦_i,z_j

Note that the wavenumber indices of the coefficients
𝖰_𝐦_i,𝐦_j,… and
𝖫_𝐦_i,… implicitly have a
MOD[𝐦_i] acting on them. Then, using the convention that sums over a basis index are implied when the associated vertical index is repeated (e.g. sum over j for repeated z_j) and that all mean-field Fourier modes
(including differences such as 𝐦_i-𝐦_j) are
explicitly written as 0, the full-basis quasilinear
equations are:

    ∂_t q (0,z_i) = 
         𝖰_0,𝐦_j,z_i,z_j,z_l q(𝐦_j,z_j)q(-𝐦_j,z_l) 
          + 𝖰_0,0,z_i,z_j,z_l q(0,z_j) q(0,z_l) 
          + 𝖫_0,z_i,z_jq(0,z_j) + 𝖢_0,z_i
       ∂_t q (𝐦_i,z_i) = 
         [gray]𝖰_𝐦_i,𝐦_j,z_i,z_j,z_l q(𝐦_j,z_j)q(𝐦_i-𝐦_j,z_l)
          + 𝖰_𝐦_i,𝐦_i,z_i,z_j,z_l q(𝐦_i,z_j) q(0,z_l) 
          + 𝖰_𝐦_i,0,z_i,z_j,z_l q(0,z_j)q(𝐦_i,z_l) 
          + 𝖫_𝐦_i,z_i,z_jq(𝐦_i,z_j) + 𝖢_𝐦_i,z_i

where the omitted nonlinear term is indicated by the strikethrough.  This expression of the quasilinear approximation in terms of Fourier modes more readily reveals a key feature and, depending on the circumstance, shortcoming of using spatial averages to define QL.  Assume that 𝖢_𝐦_i,z_i = 0, as is true for all cases studied in this paper.  Noting that all terms of the HQL eddy-field evolution (equation <ref>) are linear in q'(𝐦_n) and all terms in the HQL mean-field evolution (equation <ref>) are constant or bilinear in q'(𝐦_n), the HQL equations are seen to be invariant under an independent phase shift <cit.> for each set of all eddy modes of a specific horizontal wavenumber, 𝐦_n.  It follows that the relative phases of horizontal Fourier modes in the eddy field must only depend on the initial conditions; if these initial conditions are randomly sampled, the phases of different Fourier modes, being blind to one another, will be uncorrelated, and it is unlikely that localized structures will appear; rather, flow features in the horizontal direction (or, by extension, any homogeneous dimension of averaging used for QL) will be strictly wavelike and nonlocal, something that will manifest visually distinctly in HQL flow fields observed in the results (section <ref>).  


A reduced model may be derived by projecting the dynamical equations
of motion (equation <ref>) onto a truncated reduced
basis.  Let the POD field amplitude be written:
 
    ϕ_α = ϕ_𝐦_α,k_α = ∑_j q_𝐦_α,z_jΦ_𝐦_α,k_α,z_j^*

then the Galerkin projection of the equations of
motion onto this basis is

    ∂_t ϕ_𝐦_α,k_α   = 
       Φ_𝐦_α,k_α,z_i^*𝖰_𝐦_α,𝐦_σ,z_i,z_j,z_lΦ_𝐦_σ,k_σ,z_j
                   Φ_𝐦_α - 𝐦_σ,k_β,z_lϕ_𝐦_σ,k_σϕ_𝐦_α - 𝐦_σ,k_β
        + Φ_𝐦_α,k_α,z_i^* 𝖫_𝐦_α,z_i,z_jΦ_𝐦_α,k_σ,z_jϕ_𝐦_α,k_σ
        + Φ_𝐦_α,k_α,z_i^* 𝖢_𝐦_α,z_i

where summations are now implied when either vertical
index is repeated within a term (e.g., sum over j for repeated
k_j or z_j).  By defining the truncated evolution coefficients as

    𝖰̃_𝐦_α,𝐦_σ,k_α,k_σ,k_β  = 
           Φ_𝐦_α,k_α,z_i^* 𝖰_𝐦_α,𝐦_σ,z_i,z_j,z_lΦ_𝐦_σ,k_σ,z_jΦ_𝐦_α-𝐦_σ,k_β,z_l
       𝖫̃_𝐦_α,k_α,k_σ  = Φ_𝐦_α,k_α,z_i^* 𝖫_𝐦_α,z_i,z_jΦ_𝐦_α,k_σ,z_j
       𝖢̃_𝐦_α,k_α = Φ_𝐦_α,k_α,z_i^* 𝖢_𝐦_α,z_i

the truncated evolution equations may be succinctly expressed as

    ∂_t ϕ_0,k_α = 
         𝖰̃_0,𝐦_β,k_α,k_β,k_γϕ_𝐦_β,k_βϕ_-𝐦_β,k_γ
          + 𝖰̃_0,0,k_α,k_β,k_γϕ_0,k_βϕ_0,k_γ
          + 𝖫̃_0,k_α,k_βϕ_0,k_β + 𝖢̃_0,k_α
       ∂_t ϕ_𝐦_α,k_α = 
         [gray]𝖰̃_𝐦_α,𝐦_β,k_α,k_β,k_γϕ_𝐦_β,k_βϕ_𝐦_α-𝐦_β,k_γ
          + 𝖰̃_𝐦_α,𝐦_α,k_α,k_β,k_γϕ_𝐦_α,k_βϕ_0,k_γ
          + 𝖰̃_𝐦_α,0,k_α,k_β,k_γϕ_0,k_βϕ_𝐦_α,k_γ
          + 𝖫̃_𝐦_α,k_α,k_βϕ_𝐦_α,k_β + 𝖢̃_𝐦_α,k_α

The horizontal strikethrough indicates the terms omitted in the
reduced quasilinear approximation.








§ NUMERICAL IMPLEMENTATION



In order to compute the reduced basis for a specific, developing flow
case, an ensemble of 24 members 
was run using full-basis quasilinear (FQL) direct numerical
simulations.  Data spread evenly across all times in the flow
evolution and ensemble members is used to compute the POD modes.
Then, the energetically optimized truncated POD basis retaining
specific numbers of modes is used in a quasilinear reduced-order model
(RQL), a Galerkin Projection of the Boussinesq equations (BGP), as
described in equations <ref> and
<ref>.  For a sampling of truncated-basis sizes,
an ensemble of 12 RQL members are then compared against a reference
ensemble of 24 FQL members.

Because of the novelty of the modeling approach, this initial research
seeks a minimal and simple solver framework.  Harmonic diffusivity is
used.  Advection schemes are centered, second-order.  This caps the
maximum order of interactions at 2, simplifying quasilinearization in
all simulations.  Note that upwind advection schemes that use absolute
values of flow variables cannot be quasilinearized.



 §.§ Full-Basis Simulations


Full-basis nonlinear (FNL) and quasilinear (FQL) simulations were run
on a modified version of the MITgcm <cit.>, a finite-volume
atmospheric and oceanic fluid dynamics solver. Flow variables are
colocated in time and integrated using a second-order Adams-Bashforth
scheme. The surface boundary condition is a rigid lid.



 §.§ Reduced-Basis Determination




In order to compute the POD modes from flow data, the flow fields must
be expressed in consistent units.  This is accomplished by means of the
available potential energy (APE) in the flow

    PE =    - ∑_𝐫 g α T(𝐫) z  Δ V 
    
    APE =    PE - PE_min

where g is the acceleration of gravity at Earth's surface, α
is the coefficient of linear thermal expansion of sea water, Δ V is
the volume of each grid cell, z<0 is the fluid depth, and PE_min
is the potential energy of the temperatures sorted with lower
temperatures at greater depths. This method performed as well or
better than an alternative approach in which the temperature was
converted to a velocity such that the variance of the scaled
temperature was equal to the mean of the variances of the velocity
components. 





The POD modes are computed separately for each independent horizontal
Fourier mode using the method of snapshots (MOS) <cit.>.
Because the MITgcm's pressure solver does not yield a
divergenceless flow to numerical precision, flow snapshots are
first projected onto a divergenceless basis to entirely eliminate the
divergent modes from the decomposition.  Enough flow samples are used
to ensure at least 85% of the truncated basis subspace is
statistically converged (see figures
<ref> and
<ref>). Formally, the
convergence is computed using a convergence factor, f_c, defined as
the normalized inner product of the subspaces spanned by members of an
ensemble of bases:

    f_c = 1/N_b(N_b-1)N_r∑_αβ∑_i,j| Φ^α_i ·Φ^β_j |^2


    Φ^α_i ·Φ^β_j = δ̃_𝐦_i-𝐦_j∑_z Φ^α_i,z(Φ^β_j,z)^*

where Φ^α_i is the ith element of
basis α, N_r is the number of elements in the reduced basis,
and N_b is the number of bases in the basis ensemble. 

The cases studied here all involve developing turbulence and momentum
as well as buoyancy entrainment in the oceanic mixed-layer.  Because
the statistics are not converged, samples of flow states at  21 or 37
(depending on the case) evenly spaced times are used to determine the
bases.  Because of this, the bases retain memory of the different
states that the flow may manifest as it develops and excludes states
that should not be part of the flow solution at any time.  It is
useful to think of the basis as encoding part of the dynamics.





 §.§ Reduced-Model Implementation
 

This section serves as an introduction to the piece of software,
Boussinesq Galerkin Projection (BGP), the reduced order model which is
used in the present study and intended for use in subsequent
publications. BGP is a PDE solver that discretizes symbolic equations
of motion up to quadratic order, projects them onto a provided
truncated basis, and integrates the resulting reduced equations
explicitly forward in discrete time.  The discretization is written in
such a way that, when provided with a complete basis (on a necessarily
very small problem), BGP will be in agreement with the MITgcm up to
the precision of the latter's pressure solver. The provided basis must
be divergence-free, which allows for the reduced model to run without a pressure
solver, and defined in Fourier space, which allows for easy
quasilinearization of the quadratic term.

The computational complexity of BGP using a truncated basis with
horizontal Fourier modes (POD) may be estimated for both the model
execution and the coefficient computation (which must be computed for
every unique basis and set of flow parameters) for the RNL and RQL
solvers.  Scaling parameters will be the number of retained modes,
N_r, the number of vertical discretizations and 3D flow fields in
the underlying full-basis model, N_zf, and the number of points in
the stencil of the discretization of the quadratic term of the
underlying full-basis model, N_s. A key assumption that will be
utilized here is that the distribution function of the truncated flow
representation scales linearly with N_r such that the number of
vertical modes for a given horizontal Fourier mode scales as
N_r^1/3.  If N is the number of modes in the underlying
full-basis model, then this assumption will break down as N_r
→ N; however, this project is primarily concerned with
model reductions of N_r ≪ N.  Even for N_r ≪ N, the assumption
may still break down if the distribution function of the truncated
flow representation changes shape substantially with N_r The
following scalings are verified using runtimes from BGP in figure
<ref>:



  * the coefficient determination for RNL is 𝒪(N_r^7/3 N_zf N_s^2)

  * the coefficient determination for RQL is 𝒪(N_r^5/3 N_zf N_s^2)

  * the execution for RNL is 𝒪(N_r^7/3)

  * the execution for RQL is 𝒪(N_r^5/3)











§ TURBULENCE EXAMPLES



In order to investigate the feasibility and effectiveness of a
quasilinear reduced model for the purpose of modeling SGS
planetary turbulence, two submesoscale or smaller oceanographic flows
are used as case studies: surface-forced thermal convection and
Langmuir turbulence.  The former case is horizontally isotropic, up to
a grid-representation, while the latter is anisotropic.  Crucially,
both cases are formulated to be statistically horizontally
homogeneous, which will permit easy computation of the mean fields and
considerably simplify the POD computations and BGP execution.

The viscosity is 5000× the seawater value, and the diffusivity
is set equal to the viscosity (rather than less by about a factor of
7).  Thus, the Kolmogorov dissipation scale is roughly 600×
larger for this simulation than it would be for seawater which permits
representing the flows with a manageable grid and number of model,
i.e., an “eddy” viscosity and diffusivity closure is used.  Noting
this, it should be acknowledged that the full-resolution simulations of
these systems that are run on the MITgcm are not, strictly speaking,
direct numerical simulations of seawater.  Boundary layer problems at this resolution are
more typically managed with Large-Eddy-Simulation (LES) closures
<cit.>, but such an approach introduces
additional quadratic or higher order dependences and qualitatively
similar behavior to the eddy viscosity.  It is preferable, then, to
maintain the linear nature of the dissipation when modeling SGS
turbulent dissipation in QL dynamics. 

The FNL simulation may be considered to be a type of LES
despite the simple linear closure used here because the viscosity
and diffusivity are much larger than the seawater values and the
Prandtl number is one. LES closures generally attempt to dissipate
small-scale energy while preserving the energy cascade and other
dynamics present above the grid scale, when compared with a higher
resolution simulation of the same flow. With QL dynamics, however,
there is no energy cascade, so there can be no direct energy
backscatter; rather, all energy flows through the mean field, and as
long as the flow at the dissipation scales does not affect the
mean-field substantially, there will not be a substantial impact at
intermediate scales.  Therefore, the QL flow solution is not sensitive
to the relative amount of dissipation among different scales provided
the net amount of dissipation is approximately preserved. There is
also some flexibility in the total rate of turbulent energy
dissipation because these cases are not in a steady state and the
tendency of the horizontal mean-field is expected to balance much of
the large-scale forcing.  Roughly speaking, nonlinear LES closures
tend to dissipate energy more weakly at a given scale or at smaller
scales in comparison to harmonic viscosities and diffusivities, or to
be flow-aware selective as to where and when large effects are
needed.  Here, in the FNL case, relatively weak coefficients of
harmonic dissipation are used (gridscale Reynolds and Péclet numbers
are 𝒪(8)).  It should be stressed, however, that the
primary task here is not having an exhaustively accurate FNL
simulation, but in comparing the much larger differences between the
QL and NL runs and between the reduced-basis and full-basis
runs–particularly with QL runs. At this level of precision, the
differences between the modeling systems are much larger than
differences that result from changing subgrid schemes or parameter
values in the FNL run.



 §.§ Thermal Convection





The surface-forced-thermal-convection case (figure
<ref>) is formulated with an initial uniform
stable stratification that is cooled from the surface
<cit.>.  This causes the mixed layer to deepen.  The
nonhydrostatic Boussinesq equations that govern this case are

    ∇·𝐮   = 0 
    ∂_t 𝐮 +
    (𝐮·∇)𝐮 +
    𝐟×𝐮   = - 1/ρ_0∇
    p + 𝐠α(T - T_0) 
                    + ν∇^2 𝐮
    ∂_t T + (𝐮·∇)T
       = Q/c_p ρ_0 + κ∇^2 T

where 𝐮 = {u,v,w} is the 3-dimensional
velocity vector, T is the temperature field, t is time,
𝐟 = f 𝐳̂ is the Coriolis parameter, ρ_0
is the background density, p is the pressure, 𝐠 = g
𝐳̂ is the local acceleration of gravity, α is
the coefficient of thermal expansion, T_0 is the background
temperature, ν is the harmonic kinematic viscosity, κ is
the coefficient of thermal diffusivity, Q is the volumetric heating,
and c_p is the specific heat of sea water at constant pressure.
This case is defined on a regular, horizontally periodic Cartesian
grid with the number of cells and grid spacing per dimension
represented as {N_x,N_y,N_z} and
{Δ_x,Δ_y,Δ_z}, respectively. The domain
floor and surface have rigid, slip boundary conditions.  In the discrete model, Q = δ_k,0q/Δ_z, where k is the vertical layer index and q is
the surface heat flux.

Introducing the following non-dimensional parameters

    Ra^*    = g α q L_m^4/νκ^2 ρ_0 c_p
    
    Pr    = ν/κ

where Ra^* is the modified Rayleigh Number, appropriate for a
prescribed heat flux as opposed to a prescribed temperature at the
boundaries, Re is the Reynolds Number, Pr is the Prandtl Number,
and L_m is the mixed-layer depth.  Velocity and temperature scales
may be introduced as in <cit.>, allowing for the definition
of the Reynolds Number, a convective Froude Number, and the
conventional Rayleigh Number based on a characteristic temperature
difference, T_*:

    w_*    = √(α g q L_m/c_p ρ_0)
    
    T_*    = √(q^2/c_p^2 ρ_0^2 α g L_m)
    
    Fr    = w_*/N L_m
    
    Ra    = α g T_* L_m^3/νκ
    
    Re    = w_* L_m/ν

where N = √(α g ∂T/∂ z)
is the Brunt-Väisälä frequency, and T is the
background temperature profile.  Using these definitions, the case
studied has the following parameters:




   {N_x,N_y,N_z}       = {128,128,48} 

   {Δ_x,Δ_y,Δ_z}       = {10,10,10} m 

   f       = 0 

   ρ_0       = 1035   kg m^-3 

   g       = 9.81   m s^-2 

   α      = 2 ×10^-4   C^-1 

   T_0       = 20   C 

   ν      = 5 ×10^-2   m^2 s^-1 

   κ      = 5 ×10^-2  m^2 s^-1 

   c_p       = 3994   J kg^-1 C^-1 

   q       = -800   W m^-2 

   N^2        = 5.45 ×10^-7   s^-2 

   w_*       = 5.67 ×10^-2   m s^-1 

   T_*       = 3.41 ×10^-2   C 

   Ra^*       = 1.61 ×10^8 

   Ra       = 2.96 ×10^5 

   Pr       = 1 

   Fr       = 0.160 

   Re       = 544 
Here, N is the Brunt-Väisälä frequency.  Using these
parameters, the mixed layer depth, initially zero, will reach the
bottom of the domain at approximately t = 40 hours, at which point
the execution is stopped.  



 §.§ Langmuir Turbulence



The Langmuir-turbulence forcing is based on <cit.>.  The
domain size and resolution have been reduced in this work, while most
other case parameters are identical (see figure
<ref>).  Langmuir circulation is modeled using
the phase-averaged Craik-Leibovich equations <cit.>.  The
nonhydrostatic equations governing this case are

    ∇·𝐮 = 0 
       ∂_t 𝐮 + (𝐮·∇)𝐮 + 𝐟×( 𝐮 + 𝐮_𝐬) = 
                        - 1/ρ_0∇( p + ρ_0/2(|𝐮 + 𝐮_𝐬|^2 - |𝐮|^2 )) 
                        + ν_h ∇_h^2  𝐮 + ν_z ∇_z^2  𝐮
                        + 𝐮_𝐬×∇×𝐮 + τ/ρ_0
                        + 𝐠α(T - T_0) 
       ∂_t T + ((𝐮 + 𝐮_𝐬)·∇)T = Q/c_p ρ_0 + κ_h ∇_h^2 T + κ_z ∇_z^2 T 
       𝐮_𝐬 = u_s e^2 k_s z𝐱̂

where 𝐮_s is the Stokes drift velocity, τ
= τ𝐱̂ is the surface wind stress, k_s is the
Stokes wavenumber that sets the depth scale of the effect,
∇_h^2 are the horizontal components of the vector Laplacian,
∇_z^2 is the vertical component of the vector Laplacian,
ν_h and ν_z are the horizontal and vertical kinematic
viscosities, and κ_h and κ_z are the horizontal and
vertical thermal diffusivities.

The parameters used for this case are

   {N_x,N_y,N_z}       = {96,96,72} 

   {Δ_x,Δ_y,Δ_z}       = {3,3,0.9}   m 

   f       = 1 ×10^-4   s^-1 

   ρ_0       = 1035   kg m^-3

   g       = 9.81   m s^-2 

   α      = 2 ×10^-4   C^-1 

   T_0       = 20   C 

   ν_h       = 2 ×10^-3   m^2 s^-1 

   κ_h       = 2 ×10^-3   m^2 s^-1 

   ν_z       = 1 ×10^-3   m^2 s^-1 

   κ_z       = 1 ×10^-3   m^2 s^-1 

   c_p       = 3994   J kg^-1 C^-1 

   q       = -5   W m^-2 

   N^2        = 1.936 ×10^-5   s^-2 

   τ      = 0.037   N m^-2 

   u_s       = 0.068   m s^-1 

   k_s       = 0.105   m^-1 

   L_m        = 32.1   m 

   u_*        = 6.0 ×10^-3   m s^-1  

   Re_*        = 191.9 

   La       = 0.2965 
Here, L_m is the initial mixed-layer depth, above which the fluid
temperature is initially uniform, u_*= τ^1/2ρ_0^1/2 is the friction velocity, Re_* = u_* L_m/ν_z and the Langmuir number is defined as

    La = √(τ^1/2/ρ_0^1/2u_s)

 The surface cooling is included for
consistency with <cit.> where it is used to help
destabilize the flow; it does not have a significant impact on the
dynamics otherwise, so the Rayleigh and Prandtl number are not important.  The simulation is stopped at 72 hours.





























§ RESULTS







 §.§ Full-Basis Results



Figure <ref> shows horizontal cross sections of specific instances of the vertical velocity field of both the thermal-convection case (TC) and the Langmuir-turbulence case (LT) for FNL, FQL, and various truncations of RQL.  Figure <ref> compliments this with ensemble-averaged horizontal Fourier spectra of the flow fields taken at the same times.  The full-basis results in the top two rows of these figures will first be addressed.  

For TC vertical-velocity fields (top-left of figure <ref>), FQL has different horizontal scales and shapes than FNL.  FNL has geometric patterns of upwelling delineated by relatively intense, spoke-like plumes.  FQL, on the other hand, has very regular, horizontally wave-like plumes that are almost homogeneous in the vertical direction within the mixed layer.   The amplitude of the velocities is also observed to be substantially exaggerated in the FQL case.   Additionally, there is an apparent sign-reversal symmetry within the vertical velocity field of FQL.  This is because the Fourier modes have phases that are independent of one another, as described in section <ref>.  For a given wavenumber amplitude, there are many such independent Fourier modes on or near an annulus in wavenumber space of that radius; these must be approximately equally excited due to horizontal isotropy, creating an apparent sign-inversion symmetry within a single snapshot of the flow. 

For LT FNL vertical-velocity fields (top-right of figure <ref>), alternating bands of upwelling and downwelling are observed, broken up by small-scale turbulent features.  Like with TC, although less noticeably, the downwelling part of the structures are narrower and more concentrated.  The rolls have been rotated away from the direction of the wind and waves by the Coriolis force; the angle of rotation is dependent on depth.  FQL looks a bit like it could be a filtered version of the FNL solution (unlike TC); however, again, there is an apparent sign-inversion symmetry that is not present in FNL.  As with TC, the vertical velocity field is overly strong, likely due to a lack of small-scale turbulence acting on the primary instability.  

The Fourier spectra in the top two rows of figure <ref> reveal that in FQL, for both cases, the modes that are populated are near those which are the most energetic in FNL.  Outside of this concentrated region of wavenumbers the energy falloff is steep.  This reflects the lack of a cascade to spectrally local modes.  

The top row of figure <ref> shows vertical profiles of horizontal and ensemble-averaged mean fields for FNL and FQL, (as well as several truncations of RQL, to be discussed in section <ref><ref>).  FQL is able to reproduce the FNL fairly well.  The horizontal mean fields are particularly important in the quasilinear approximation as they are the only way turbulent energy can move among horizontal Fourier modes, so it is expected that the rough agreement between FQL and FNL mean profiles is an important foundation for making further approximations in the form of model reduction. 

In contrast to the mean-field profiles (figure <ref>), the second-order FQL profiles (figures <ref>, <ref>, and <ref>) generally do not exhibit good agreement with the FNL profiles.  For the thermal-convection case, the only nonlinear profile that is well represented by FQL is the vertical turbulent temperature transport, w'T', figure <ref>.  This is likely because the transport of temperature is what is driving the flow; a mean-field temperature anomaly builds at the surface due to the prescribed heat flux and is preferentially transported downwards due to the eddy field.  This term dominates the turbulent-kinetic-energy production, which is also accurately represented although not depicted.  Note that other off-diagonal Reynolds stresses are exactly zero in the limit of a large ensemble average because of the horizontal isotropy of this case, so there are relatively few non-trivial second-order vertical statistical profiles to observe in TC.  For LT, none of the second-order FQL profiles are well-represented.





 §.§ Reduced-Model Fields and Spectra



Reduced-basis fields and horizontal Fourier spectra for TC and LT are provided in figures <ref> and <ref>, respectively, along with those of the full bases.    At a glance, the behavior of TC and LT are substantially different under the model reduction. In the former case, modes are populated at slightly smaller wavenumbers than FQL while in LT at truncations larger than 50, energy tends to be in much larger wavenumbers than FQL.  The RQL flow solution and statistics must eventually converge upon those of FQL as the basis size is further increased.  The convergence of the Fourier energy footprint is therefore nonuniform, where smaller basis retentions can actually look more like FQL than the larger basis retentions.  It will turn out to be the case that this nonuniform convergence can be advantageous in that it allows for better representation of the FNL statistical profiles.   The nonuniform convergence of the RQL flow statistics on those of FQL is explored more thoroughly in a small test case in section <ref><ref>. 

In TC RQL 50, the observed rows of plumes are due to an insufficient number of retained eddy-modes on the dominant wavenumber (see bottom left of figure <ref>).  The basis was selected in a manner that observes discrete horizontal isotropy exactly, so the non-zero modes show an eight-fold symmetry; however, with fifty modes, the truncation is so severe that the annulus is not complete and modes along the axes are omitted.  Because so few modes are excitable near this wavenumber, their amplitudes do not average out to form smooth structures that reflect the symmetry of the problem, rather, one or two may randomly dominate at a given time, giving rise in some instances to rows, in others to regular triangular patterns.  Despite the qualitative shortcomings, the vertical profile resembles the FQL solution, exaggerating the plume depth slightly.  In LT RQL 50, the breaking of horizontal isotropy allows for 50 modes to better capture the essential shape of the windrows than with TC; however, the velocity field is too strong and penetrates through the stratified layer.  

At 500 modes, both TC and LT show improvement in the vertical cross sections of the flow fields (figure <ref>), a reflection of the larger number of accessible modes, as seen in the 2D Fourier spectrum in figure <ref>.  The horizontal cross sections of both cases reflect RQL's shift to different horizontal scales compared with FQL; in TC RQL 500, the solution has similar horizontal structure as FQL except with slightly larger length scales, while in LT RQL 500, the windrows are much more narrow.   This trend is more pronounced at 5000 modes.  The displacement of the excitation of horizontal Fourier modes at intermediate basis truncations in RQL is explored more thoroughly in section <ref><ref>. 





 §.§ Reduced-Model Statistical Profiles



The mean-field profiles are depicted in figure <ref> along with their representation error:

    e_i = (<q_i, FNL> - <q_i, RQL>) ·<q_i, FNL>/<q_i, FNL> ·<q_i, FNL>

Here, the mean field is <q(z,t)>, <> indicates an ensemble average, the subscript runs over the fields {u,v,w,T}, t runs over four evenly sampled times throughout the flow evolution, z runs over the vertical position, and the dot product is defined over both z and t.  The top-rows of figures <ref>, <ref>, and <ref> depict vertical profiles of horizontal and ensemble averaged turbulent transports and energy terms, while the bottom rows show the RQL model's representation error, e_ij of the FNL statistical profiles, computed as a normalized inner product of the Reynolds stresses <q_i q_j> = τ_ij = τ_ij(z,t):

    e_ij = (τ_ij, FNL - τ_ij, RQL) ·τ_ij, FNL/τ_ij, FNL·τ_ij, FNL

Note that errors less than 0.3 usually indicate a good qualitative representation of the vertical profile. 

At small basis sizes (200 and 1000), the mean-field thermal-convection temperature profile, depicted on the left of figure <ref>, as well as the vertical temperature transport, w'T' (right of figure <ref>), are represented well by RQL.  At 5000 modes, however, unrealistically strong and deep profiles of w'w' and w'T' in figures <ref> and <ref>, respectively, have mixed the domain entirely by the time of 20 hours, demonstrating a deviation from uniform convergence towards the FQL statistics.  An explanation and solution for this anomalous mixing using a modified POD basis is given in section <ref><ref>.

Unlike the thermal-convection case, the Langmuir-case momentum mean field profiles in figure <ref> do not exhibit a trend of worse representation at larger basis sizes. The FQL v-velocity profile appears qualitatively consistent with the FNL profiles; however, the u-velocity component exhibits errors on the order of the magnitude of the FNL mean field (meaning e_u ≈ 1). Note that for none of the mean-field profiles does RQL outperform FQL at reproducing the FNL mean fields.  This may be unsurprising given the degree of reduction of the model size; however, it will turn out that RQL will sometimes perform better than FQL for second-order profiles. 

The Langmuir-turbulence case, on the other hand, exhibits good representation of several turbulent transport and energy profiles with RQL; however, it requires a larger portion of retained basis modes to do this, about 0.2% of modes, or about a factor of 10 more than the thermal-convection case. Several of the RQL turbulent transport and energy profiles (specifically u'u', w'w', u'v' and v'w'), appear to initially converge on the FNL profiles (as opposed to FQL), something which was not observed in the mean-field profiles in figure <ref>.  The potential of RQL to reproduce second-order FNL statistics would be a huge boon to sub-grid-scale modeling within the ocean-surface boundary layer, so the mechanism for this convergence is scrutinized in section <ref><ref>.  It can also be noted that the relatively good performance of RQL at reproducing second-order FNL profiles compared with first-order mean profiles indicates the mean fields may not be as essential as would have been expected given their important role as an energy pathway in the quasilinear approximation.





 §.§ Reduced Basis Excitation Displacement



In order to understand the convergence of RQL on FNL statistics in LT as well as the degradation of performance of RQL TC at intermediate basis truncations, it is helpful to first identify the cause and impact of the change in excited horizontal Fourier wavenumbers that occurs in RQL relative to FQL (see figure <ref>).  Recall that relative to FQL, RQL 5000 prefers smaller wavenumbers in TC and quite a bit larger in LT. To probe the mechanism for this spectral excitation displacement, simulations were run to determine the growth rate of quasilinear instabilities around a fixed FQL mean field for different basis truncations. Based on these numerical experiments, the peak quasilinear instability is displaced to different horizontal scales because the truncation of these excited horizontal modes prevents the redistribution of energy to different vertical profiles within that horizontal mode. This then prevents the draining of energy that would otherwise occur in the full-basis solution. In the LT fixed-mean-field experiments, high-wavenumber modes are driven by the linear Stokes instability. In the full-basis flow solution, the vertical profile of this mode evolves through an interaction with the mean field, and backscatter and vertical viscosity both act to dissipate this energy such that the mode never becomes substantially energized. In the 5000-mode version, however, the representation of these horizontal Fourier modes is limited and energy is not able to vertically redistribute in order to drain through these mechanisms. Similarly, in TC the smaller-wavenumber modes are excited via scattering off of the mean field and efficiently dissipated via backscatter, viscosity, and thermal diffusion in the full-basis solution, with the truncated basis (1000+) suppressing this pathway.  

Having identified the cause of these shifted peak excitations to be an insufficient representation of the associated horizontal Fourier modes, next consider their impact on the mean-field and vertical structure of the flow.  Recall that the POD bases are computed using snapshots of FQL flow solutions.  In FQL, boundary layer structures are simple, spectrally focused, and uniform throughout the mixed layer (relative to FNL).  The horizontal dimension of these flow structures scales monotonically with the mixed-layer depth.  When the basis truncation of a horizontal Fourier mode is very limited, such as at the horizontal wavenumbers of the displaced peak excitations in LT and TC, the included vertical modes will allow for only a small variation of mixed-layer depths around the preferred depth of the original FQL solution.  Because of this, in LT the displaced peak-excitation Fourier modes of a high-wavenumber create structures of a shallower profile than the overall MLD (which is controlled by some of the less energetic and smaller wavenumber modes).  The interior of the mixed layer becomes structurally less homogeneous, and therefore more FNL-like (bottom-right of figure <ref>).  For TC, the displaced peak-excitation horizontal Fourier modes artificially push the mixed-layer beyond its correct size, causing the RQL profiles to be a much worse representation of the FQL results (Fig. <ref>).  

The reason that TC and not LT exhibits displaced mode excitation at larger scales is because in TC the representation of small-wavenumber modes gradually trails off at small wavenumbers in the FQL solution while for LT they are continuously represented up to a sharp cutoff, missing only a single discrete wavenumber next to the mean field (see figure <ref>.)  A means of addressing the inaccuracies associated with large-scale displaced RQL mode excitation while preserving the benefits of the small-scale modes is to allow the reduced-basis solver to use only those partially represented modes that are approximately of equal or larger horizontal wavenumber than the peak excitations. This basis modification could be efficient to implement because it can take effect at runtime, which stands in contrast to modifications to the POD formulation, such as different energy scalings of the temperature, which must be configured in advance of basis determination.

To test this idea, the reduced basis was further truncated at poorly represented wavenumbers that are smaller than the modes of peak representation (see figure <ref>). Poorly represented is taken as 85% or less of the peak represented mode. This was done with the understanding that a dynamically adaptive solver would perform as well or better than this test with similar computational expense. The results of this test, shown in figures <ref> and <ref>, show that the technique can avoid the most severe inaccuracies of the truncated large-scale horizontal modes and also preserve the approximately accurate evolution of the deepening mixed layer.  






 §.§ Miniature Langmuir-Case Convergence Study



In order to examine more closely the convergence properties of RQL on FQL, a miniature test system containing 10 × 10 × 10 grid cells was simulated that could be run using the Galerkin Projection of the equations of motion all the way up to complete resolution of the underlying basis.  The system parameters are exactly the same as the Langmuir turbulence case (see section <ref>), with the exception of the mixed-layer depth, which is initialized at 5 meters.  The grid spacing is the same in both cases, meaning the domain only extends to 9 meters depth instead of 72 meters depth for the larger case.  The resulting flow field is far from realistic and the mixed layer quickly reaches the bottom of the domain; however, the nonlinear and quasilinear flow solutions do have substantially different vertical statistical profiles prior to that, during the first 3 hours of the flow evolution.  384 ensemble members were run of FNL, FQL, and RQL; the simulations were run until 10 hours into the flow evolution. 

The miniature test system indicates that for some turbulent transport profiles, the RQL model trends towards the FNL profiles at intermediate basis sizes, for instance, the vertical turbulent x-component velocity transport u'w' (left side of figure <ref>).  This quantity is responsible for transporting much of the turbulent energy in the flow introduced near the surface to greater depths and Monin-obukhov similarity theory constrains this quantity heavily in this case.  Perhaps it is for these reasons that u'w' is represented relatively well in both the miniature test case and in the larger Langmuir-turbulence case. 

Other RQL turbulent profiles are representative of the FNL profiles at intermediate basis retention because the former incidentally sweeps through the FNL profile as with w'w' in the center column of figure <ref>.  Finally, some profiles do not exhibit agreement between the FNL and RQL at intermediate basis retention, as in v'w' in the right column where the RQL reproduces only part of the FNL profile at a given truncation size. The magnitude of this transport is considerably smaller than u'w', perhaps suggesting why it might not be accurately represented.  Many of these turbulent profiles do not exhibit uniform convergence towards the FQL versions.  This can be seen in v'w' in the right column of figure <ref> most clearly among the three statistical quantities depicted.

Overall, the miniature Langmuir-turbulence test case demonstrates that the intermediate convergence exhibited by many of the statistical profiles in the larger Langmuir case towards the FNL profile (see figures <ref> and <ref>) as well as the nonuniform convergence found in all statistical profiles in the thermal-convection case in figure <ref> occurs in other cases and is not an unusual feature or error.











§ DISCUSSION





One takeaway from this work is that some turbulent-transport and
energy profiles are well represented in the reduced models, while
others are not.  The ones that are well represented seem to be those
which are most important for transporting energy from externally
forced boundary conditions into the flow interior.  It would be
useful to be able to predict which of these transports would
be most accurately represented without necessarily referring back to
an ensemble of full-basis simulations.  Monin-Obukhov similarity
theory or <cit.> budgets might be revealing, but more
simulations with both convective and stress forcing are needed to
evaluate the robustness of this hypothesis.

As an idea for what might indicate the importance of second-order
statistical transports, a 4-dimensional Reynolds-stress tensor was
computed using the temperature scaling described in section
<ref><ref> as well as the APE definition in
equation <ref> and normalized using the variances of all the
elements at different depths and times:

    ϵ_ij = ∑_z,tq_i'q_j'(z,t)^2/∑_{ij}∑_z,tq_i'q_j'(z,t)^2
 
where {ij} is the set of pairs of velocity-temperature
components in consideration (energy and transport elements are treated
separately for best results, so this will be either i=j or i j).
This was then plotted against the representation error, e_ij (see
equation <ref>), for a variety of different cases.  The
results, shown in figure <ref>, indicate a weakly
negative correlation between the error and the scaled Reynolds stress
amplitude.  Perhaps a more reasoned approach to determining the
“importance” of these statistics might reveal a stronger
correlation.  











§ CONCLUSION



Two cases of developing ocean-surface boundary-layer turbulence were used to study the quasilinear
reduced-model representation of flow statistics.  In concrete terms,
the reduced quasilinear model represented certain second-order
turbulent transport and energy profiles in these cases qualitatively
well (less than 30% error from the full-basis nonlinear profiles)
using 𝒪(10^-3) retained POD modes.  The
profiles that were represented well appear to be correlated with their role
driving the flow dynamics.  

The case of surface-forced thermal convection in a deepening boundary
layer exhibited nonuniform convergence of the statistical
profiles towards complete representation, causing substantial inaccuracies 
in flow statistics.  This was corrected for by restricting basis representation 
at large scales in a way that can be done dynamically at runtime.  The reduced-basis
quasilinear vertical statistics profiles in the Langmuir-turbulence
case, of which there are many more than in the thermal-convection case
that are non-trivial, exhibited intermediate convergence towards the
full-basis nonlinear profiles before converging on the full-basis
quasilinear profiles.  Both of these properties were observed more
thoroughly in an unrealistic miniature test case of the Langmuir
turbulence. 

What this means in broad terms is that a QL POD reduced model performs
better than would be expected based naively on the isolated
performances of the QL approximation and POD-based reduced-modeling:
not only does QL allow for a smaller problem size and better runtime
scalings than nonlinear POD reduced model; truncation of the QL
approximation can add multi-scale structure missing from the full-basis QL
solution, often producing flow features that resemble those of
full-basis NL more so than full-basis QL. By better understanding
these synergistic effects, it may be possible to craft QL POD reduced
models into a cost-effective SGS model. 

In addition to the performance of the cases studied, the potential for
generalizing this framework can be acknowledged.  Although only two
examples of turbulence were tested of the many that must be modeled in
an ESM, it should be noted that these were the only examples studied.
Because there were no assumptions made that are specific to these
cases in developing this QL POD reduced-model framework, it is
reasonable to expect the framework to work similarly well with other
types of boundary-layer turbulence in general, potentially addressing the drawbacks
of the current approach of modeling different turbulent phenomena
separately.




§ ACKNOWLEDGEMENTS

We thank Altan Allawala, Greg Chini, and Steve Tobias for helpful conversations.  The authors are grateful to the 
Kavli Institute for Theoretical Physics for supporting the Program “Planetary Boundary Layers in Atmospheres, Oceans, and Ice on Earth and Moons” in Spring 2018 where some of this work was carried out.  The research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958, NSF ONR N00014-17-1-2963, NSF OCE-1350795 as well as by the Institute at Brown for Environment and Society.









§ PROPER ORTHOGONAL DECOMPOSITION IN HOMOGENEOUS DIMENSIONS



Begin with a field, ψ(𝐱), defined over n
continuous and statistically homogeneous periodic dimensions that span
0 to L_i and define the volume V = ∏_i^n
  L_i.  For simplicity, inhomogeneous dimensions are not written
explicitly but all variables can be a defined over any number of
inhomogeneous dimensions without loss of generality. Because the
statistics are homogeneous in 𝐱, the covariance, c_2,
indicating the second cumulant, must depend only on the field spacing:

    c_2(𝐱_1 - 𝐱_2) = ψ(𝐱_1)ψ(𝐱_2)

where the overline now specifically refers to a spatial average over
the n homogeneous dimensions.  Introducing 𝐱̅ =
1/2(𝐱_1 + 𝐱_2) and δ𝐱 = 𝐱_1 - 𝐱_2, then this average can be
written

    c_2(𝐱_1 - 𝐱_2) = 1/V∫_𝒱_x̅ d𝐱̅ψ(𝐱̅ + δ𝐱/2)ψ(𝐱̅ - δ𝐱/2)

Where 𝒱_x̅ spans the entire domain in position
space, [0,L_i) for the ith homogeneous dimension.
  Then, writing the fields as a sum of their Fourier modes

    ψ(𝐱_1)    = ∫_𝒱_k_1d𝐤_1/(2 π)^n e^i 𝐤_1 ·𝐱_1ψ(𝐤_1) 
    ψ(𝐱_2)    = ∫_𝒱_k_2d𝐤_2/(2 π)^n e^i 𝐤_2 ·𝐱_2ψ(𝐤_2)

where 𝒱_k spans the entire volume of wavenumber space,
(-∞,+∞) in each homogeneous dimension.
Plugging this into the spatial average

    c_2(𝐱_1 - 𝐱_2) = 
           1/V∫_𝒱_k_1,𝒱_k_2d𝐤_1/(2 π)^nd𝐤_2/(2 π)^n∫_𝒱_x̅ d𝐱̅
                e^i 𝐤_1 ·(𝐱̅+δ𝐱/2) e^i 𝐤_2 ·(𝐱̅-δ𝐱/2)ψ(𝐤_1)ψ(𝐤_2) 
        c_2(𝐱_1 - 𝐱_2) = 
           ∫_𝒱_kd𝐤_1 d𝐤_2/(2 π)^nδ(𝐤_1+𝐤_2) 
                e^i 𝐤_1 ·δ𝐱/2 e^- i 𝐤_2 ·δ𝐱/2ψ(𝐤_1)ψ(𝐤_2) 
        c_2(𝐱_1 - 𝐱_2) = 
           ∫_𝒱_kd𝐤_2/(2 π)^n e^i 𝐤_1 ·δ𝐱ψ(𝐤_1)ψ(𝐤_1)^*

so the covariance matrix is diagonal in Fourier space, and therefore
the POD modes are also Fourier modes in homogeneous dimensions. 









ametsoc2014




59




urlstyle

[Adcroft et al.(2018)]mitgcm
Adcroft, A., and Coauthors, 2018: Mitgcm user manual. MIT department of EAPS,
  <http://mitgcm.org/download/manual/manual_20180114.pdf>.

[Allawala and Marston(2016)Allawala, and Marston]Allawala:2016hx
Allawala, A., and J. B. Marston, 2016: Statistics of the stochastically-forced
  Lorenz attractor by the Fokker-Planck equation and cumulant expansions.
  Physical review. E, Statistical physics, plasmas, fluids, and related
  interdisciplinary topics, 94 (5), 052 218.

[Allawala et al.(2017)Allawala, Tobias,, and Marston]allawala17
Allawala, A., S. M. Tobias, and J. B. Marston, 2017: Dimensional reduction of
  direct statistical simulation. arXiv preprint arXiv:1708.07805.

[Bachman et al.(2017)Bachman, Fox-Kemper, Taylor,, and
  Thomas]bachman2017parameterization
Bachman, S., B. Fox-Kemper, J. Taylor, and L. Thomas, 2017: Parameterization of
  frontal symmetric instabilities. i: Theory for resolved fronts. Ocean
  Modelling, 109, 72–95.

[Bartello and Holloway(1991)Bartello, and
  Holloway]bartelloholloway1991
Bartello, P., and G. Holloway, 1991: Passive scalar transport in β-plane
  turbulence. Journal of Fluid Mechanics, 223, 521–536,
  10.1017/S0022112091001532,
  <http://journals.cambridge.org/article_S0022112091001532>.

[Batchelor(1952)]batchelor1952diffusion
Batchelor, G., 1952: Diffusion in a field of homogeneous turbulence: Ii. the
  relative motion of particles. Mathematical Proceedings of the
  Cambridge Philosophical Society, Cambridge University Press, Vol. 48,
  345–362.

[Bauer et al.(2015)Bauer, Thorpe,, and Brunet]Bauer_2015
Bauer, P., A. Thorpe, and G. Brunet, 2015: The quiet revolution of numerical
  weather prediction. Nature, 525 (7567), 47–55,
  10.1038/nature14956,
  <http://dx.doi.org/10.1038/nature14956>.

[Belcher et al.(2012)]belcher2012global
Belcher, S. E., and Coauthors, 2012: A global perspective on langmuir
  turbulence in the ocean surface boundary layer. Geophysical Research
  Letters, 39 (18).

[Bony et al.(2015)]Bony:2015bd
Bony, S., and Coauthors, 2015: Clouds, circulation and climate sensitivity.
  Nature Geoscience, 8 (4), 261–268.

[Bouchet et al.(2013)Bouchet, Nardini,, and
  Tangarife]bouchetnardinietal2013
Bouchet, F., C. Nardini, and T. Tangarife, 2013: Kinetic Theory of Jet
  Dynamics in the Stochastic Barotropic and 2D Navier-Stokes Equations.
  Journal of Statistical Physics, 153 (4), 572–625,
  10.1007/s10955-013-0828-3,
  <http://dx.doi.org/10.1007/s10955-013-0828-3>.

[Canuto and Minotti(2001)Canuto, and
  Minotti]canutominotti2001
Canuto, V. M., and F. Minotti, 2001: Mixing and transport in stars - I.
  Formalism: momentum, heat and mean molecular weight. Mon. Not. Roy.
  Ast. Soc., 328, 829–838.

[Chaalal et al.(2016)Chaalal, Schneider, Meyer,, and
  Marston]Chaalal:2016jx
Chaalal, F. A., T. Schneider, B. Meyer, and B. Marston, 2016: Cumulant
  expansions for atmospheric flows. New Journal of Physics,
  18 (2), 1–24.

[Craik(1977)]craik77
Craik, A. D. D., 1977: The generation of langmuir circulations by an
  instability mechanism. Journal of Fluid Mechanics, 81 (2),
  209–223.

[Doney et al.(2004)]Doney_2004
Doney, S. C., and Coauthors, 2004: Evaluating global ocean carbon models: The
  importance of realistic physics. Global Biogeochem. Cycles,
  18 (3), n/a–n/a, 10.1029/2003gb002150,
  <http://dx.doi.org/10.1029/2003gb002150>.

[Farrell and Ioannou(2007)Farrell, and Ioannou]farrell07
Farrell, B. F., and P. J. Ioannou, 2007: Structure and spacing of jets in
  barotropic turbulence. Journal of the Atmospheric Sciences,
  64 (10), 3652–3665.

[Foken(2006)]foken200650
Foken, T., 2006: 50 years of the monin–obukhov similarity theory.
  Boundary-Layer Meteorology, 119 (3), 431–447.

[Fox-Kemper et al.(2014)Fox-Kemper, Backman, Pearson,, and
  Reckinger]fox-kemper14
Fox-Kemper, B., S. D. Backman, B. Pearson, and S. Reckinger, 2014: Principles
  and advances in subgrid modelling for eddy-rich simulations. Clivar
  Exchanges, 19 (2), 42–46.

[Fox-Kemper and Ferrari(2008)Fox-Kemper, and Ferrari]fox-kemper08
Fox-Kemper, B., and R. Ferrari, 2008: Parameterization of mixed layer eddies.
  part i: Theory and diagnosis. Journal of Physical Oceanography,
  38 (6), 1145–1165.

[Hamlington et al.(2014)Hamlington, Roekel, Fox-Kemper, Julien,, and
  Chini]Hamlington_2014
Hamlington, P. E., L. P. V. Roekel, B. Fox-Kemper, K. Julien, and G. P. Chini,
  2014: Langmuir–Submesoscale Interactions: Descriptive Analysis
  of Multiscale Frontal Spindown Simulations. J. Phys. Oceanogr.,
  44 (9), 2249–2272, 10.1175/jpo-d-13-0139.1,
  <http://dx.doi.org/10.1175/JPO-D-13-0139.1>.

[Harcourt(2013)]harcourt2013second
Harcourt, R. R., 2013: A second-moment closure model of langmuir turbulence.
  Journal of Physical Oceanography, 43 (4), 673–697.

[Herring(1963)]herring63
Herring, J. R., 1963: Investigation of problems in thermal convection.
  Journal of Atmospheric Sciences, 20, 325–338.

[Holmes et al.(2012)Holmes, Lumley, Berkooz,, and Rowley]holmes12
Holmes, P., J. L. Lumley, G. Berkooz, and C. W. Rowley, 2012:
  Turbulence, Coherent Structures, Dynamical Systems and Symmetry. 2nd
  ed., Cambridge University Press, Cambridge.

[Kang et al.(2009)Kang, Frierson,, and Held]kang2009tropical
Kang, S. M., D. M. Frierson, and I. M. Held, 2009: The tropical response to
  extratropical thermal forcing in an idealized gcm: The importance of
  radiative feedbacks and convective parameterization. Journal of the
  atmospheric sciences, 66 (9), 2812–2827.

[Kolmogorov(1941)]kolmogorov1941dissipation
Kolmogorov, A. N., 1941: Dissipation of energy in locally isotropic turbulence.
  Dokl. Akad. Nauk SSSR, Vol. 32, 16–18.

[Kraus and Turner(1967)Kraus, and Turner]KrausTurner67
Kraus, E., and J. Turner, 1967: A one-dimensional model of the seasonal
  thermocline. II: The general theory and its consequences.
  Tellus, 19, 98–106.

[Krause and Raedler(1980)Krause, and
  Raedler]krauseraedler1980
Krause, F., and K. H. Raedler, 1980: Mean-field
  magnetohydrodynamics and dynamo theory. Pergamon Press.

[Large et al.(1997)Large, Danabasoglu, Doney,, and
  McWilliams]Large_1997
Large, W. G., G. Danabasoglu, S. C. Doney, and J. C. McWilliams, 1997:
  Sensitivity to Surface Forcing and Boundary Layer Mixing in a Global Ocean
  Model: Annual-Mean Climatology. J. Phys. Oceanogr.,
  27 (11), 2418–2447,
  10.1175/1520-0485(1997)027<2418:stsfab>2.0.co;2,
  <http://dx.doi.org/10.1175/1520-0485(1997)027<2418:stsfab>2.0.co;2>.

[Large et al.(1994)Large, McWilliams,, and Doney]large94
Large, W. G., J. C. McWilliams, and S. C. Doney, 1994: Oceanic vertical mixing:
  A review and a model with a nonlocal boundary layer parameterization.
  Reviews of Geophysics, 32 (4), 363 – 403.

[Leibovich(1983)]Leibovich_1983
Leibovich, S., 1983: The form and Dynamics of Langmuir Circulations.
  Annu. Rev. Fluid Mech., 15 (1), 391–427,
  10.1146/annurev.fl.15.010183.002135,
  <http://dx.doi.org/10.1146/annurev.fl.15.010183.002135>.

[Li et al.(2019)]LiReichl19
Li, Q., and Coauthors, 2019: Comparing ocean boundary vertical mixing schemes
  including Langmuir turbulence. Journal of Advances in Modeling
  Earth Systems (JAMES).

[Malkus(1954)]MALKUS:1954dh
Malkus, W., 1954: The Heat Transport and Spectrum of Thermal Turbulence.
  Proceedings of the Royal Society of London. Series A, Mathematical
  and Physical Sciences, 225 (1161), 196–212.

[Maltrud and Vallis(1991)Maltrud, and Vallis]maltrudvallis1991
Maltrud, M. E., and G. K. Vallis, 1991: Energy spectra and coherent structures
  in forced two-dimensional and beta-plane turbulence. Journal of
  Fluid Mechanics, 228, 321–342.

[Marston et al.(2016)Marston, Chini,, and Tobias]marston16
Marston, J., G. Chini, and S. Tobias, 2016: Generalized quasilinear
  approximation: Application to zonal jets. Physical Review Letters,
  116 (21), 214 501.

[Marston(2010)]marston10
Marston, J. B., 2010: Statistics of the general circulation from cumulant
  expansions. Chaos, 20 (4), 041 107.

[Marston et al.(2008)Marston, Conover,, and
  Schneider]marstonconoveretal2008
Marston, J. B., E. Conover, and T. Schneider, 2008: Statistics of an
  Unstable Barotropic Jet from a Cumulant Expansion. Journal of
  Atmospheric Sciences, 65, 1955, 10.1175/2007JAS2510.1,
  0705.0011.

[Marston et al.(2019)Marston, Qi,, and Tobias]marston14
Marston, J. B., W. Qi, and S. M. Tobias, 2019: Direct Statistical
  Simulation of a Jet in Zonal Jets: Phenomenology, Genesis, and
  Physics (arXiv:1412.0381). Cambridge University Press.

[McWilliams and Sullivan(2000)McWilliams, and
  Sullivan]McWilliams_2000
McWilliams, J. C., and P. P. Sullivan, 2000: Vertical Mixing by Langmuir
  Circulations. Spill Science & Technology Bulletin,
  6 (3-4), 225–237, 10.1016/s1353-2561(01)00041-x,
  <http://dx.doi.org/10.1016/s1353-2561(01)00041-x>.

[McWilliams et al.(1997)McWilliams, Sullivan,, and
  Moeng]mcwilliams97
McWilliams, J. C., P. P. Sullivan, and C. Moeng, 1997: Langmuir turbulence in
  the ocean. Journal of Fluid Mechanics, 334, 1–30.

[Mellor and Yamada(1974)Mellor, and Yamada]mellor74
Mellor, G. L., and T. Yamada, 1974: A hierarchy of turbulence closure models
  for planetary boundary layers. Journal of the Atmospheric Sciences,
  31, 1791–1806.

[O'Gorman and Schneider(2007)O'Gorman, and Schneider]ogorman07
O'Gorman, P. A., and T. Schneider, 2007: Recovery of atmospheric flow
  statistics in a general circulation model. Geophysical Research
  Letters, 34, 524 – 535.

[Pausch et al.(2019)Pausch, Yang, Hwang,, and Eckhardt]pausch2019
Pausch, M., Q. Yang, Y. Hwang, and B. Eckhardt, 2019: Quasilinear approximation
  for exact coherent states in parallel shear flows. Fluid Dynamics
  Research, 51, 1–15.

[Rüdiger(1989)]ruediger1988
Rüdiger, G., 1989: Differential rotation and stellar convection:
  Sun and solar-type stars, Vol. 5. Gordon & Breach, New York.

[Salmon(1998)]salmon98
Salmon, R., 1998: Lectures on Geophysical Fluid Dynamics. Oxford
  University Press, New York.

[Schmidt and Schumann(1989)Schmidt, and Schumann]schmidt89
Schmidt, H., and U. Schumann, 1989: Coherent structure of the convective
  boundary layer derived from large-eddy simulation. Journal of Fluid
  Mechanics, 200, 511–562.

[Sirovich(1989)]sirovich89
Sirovich, L., 1989: Chaotic dynamics of coherent structures. Physica D:
  Nonlinear Phenomena, 37 (1-3), 126–145.

[Skitka(2019)]skitka19_thesis
Skitka, J., 2019: Quasilinear modeling of planetary boundary-layer turbulence
  [PhD Thesis]. Providence, RI: Brown University.

[Spiegel(1962)]SPIEGEL:1962vf
Spiegel, E. A., 1962: Thermal Turbulence at Very Small Prandtl Number.
  Journal of Geophysical Research, 67 (8), 3063–&.

[Squire and Bhattacharjee(2015)Squire, and
  Bhattacharjee]Squire:2015fk
Squire, J., and A. Bhattacharjee, 2015: Statistical Simulation of the
  Magnetorotational Dynamo. Physical Review Letters,
  114 (8), 085 002–5.

[Srinivasan and Young(2012a)Srinivasan, and
  Young]srinivasan12
Srinivasan, K., and W. Young, 2012a: Zonostrophic instability.
  Journal of the atmospheric sciences, 69 (5), 1633–1656.

[Srinivasan and Young(2012b)Srinivasan, and
  Young]srinivasanyoung2012
Srinivasan, K., and W. R. Young, 2012b: Zonostrophic
  instability. J. Atmos. Sci., 69, 1633–1656.

[Stephens(2005)]2005JClimate
Stephens, G. L., 2005: Cloud Feedbacks in the Climate System: A Critical
  Review. Journal of Climate, 18 (2), 237–273.

[Stocker et al.(2013)]IPCC13
Stocker, T. F., and Coauthors, 2013: Climate change 2013: The physical science
  basis. Tech. rep. 1535 pp.

[Sullivan et al.(1994)Sullivan, McWilliams,, and Moeng]sullivan94
Sullivan, P. P., J. C. McWilliams, and C.-H. Moeng, 1994: A subgrid-scale model
  for large-eddy simulation of planetary boundary-layer flows.
  Boundary-Layer Meteorology, 71 (3), 247–276.

[Thomas et al.(2015)Thomas, Farrell, Ioannou,, and
  Gayme]thomasetal2015
Thomas, V. L., B. F. Farrell, P. J. Ioannou, and D. F. Gayme, 2015: A
  minimal model of self-sustaining turbulence. Physics of Fluids,
  27 (10), 105 104, 10.1063/1.4931776.

[Thorpe(2004)]Thorpe_2004
Thorpe, S., 2004: Langmuir Circulation. Annual Review of Fluid
  Mechanics, 36 (1), 55–79,
  10.1146/annurev.fluid.36.052203.071431,
  <http://dx.doi.org/10.1146/annurev.fluid.36.052203.071431>.

[Tobias et al.(2011)Tobias, Dagon,, and Marston]tobias11
Tobias, S. M., K. Dagon, and J. B. Marston, 2011: Astrophysical fluid dynamics
  via direct statistical simulation. The Astrophysical Journal,
  727, 127.

[Tobias and Marston(2013)Tobias, and Marston]tobias13
Tobias, S. M., and J. B. Marston, 2013: Direct Statistical Simulation of
  Out-of-Equilibrium Jets. Physical Review Letters,
  110 (10), 104 502.

[Van Roekel et al.(2012)Van Roekel, Fox-Kemper, Sullivan, Hamlington,,
  and Haney]van2012
Van Roekel, L., B. Fox-Kemper, P. Sullivan, P. Hamlington, and S. Haney, 2012:
  The form and orientation of langmuir cells for misaligned winds and waves.
  Journal of Geophysical Research: Oceans, 117 (C5).

[Vedenov et al.(1961)Vedenov, Velikhov,, and Sagdeev]Vedenov:1961us
Vedenov, A. A., E. P. Velikhov, and R. Z. Sagdeev, 1961: Quasilinear theory of
  plasma oscillations. Proceedings of IAEA Conference on Plasma
  Physics and Controlled Nuclear Fusion Research, 465–475.

[Wang et al.(2012)Wang, Akhtar, Borggaard,, and Iliescu]wang12
Wang, Z., I. Akhtar, J. Borggaard, and T. Iliescu, 2012: Proper orthogonal
  decomposition closure models for turbulent flows: a numerical comparison.
  Computer Methods in Applied Mechanics and Engineering, 237,
  10–26.








