Quantifying mixing and available potential energy in vertically periodic simulations of stratified flows

Turbulent mixing exerts a significant influence on many physical processes in the ocean. In a stably stratified Boussinesq fluid, this irreversible mixing describes the conversion of available potential energy (APE) to background potential energy (BPE). In some settings the APE framework is difficult to apply and approximate measures are used to estimate irreversible mixing. For example, numerical simulations of stratified turbulence often use triply periodic domains to increase computational efficiency. In this setup however, BPE is not uniquely defined and the method of Winters et al. (1995, J. Fluid Mech., 289) cannot be directly applied to calculate the APE. We propose a new technique to calculate APE in periodic domains with a mean stratification. By defining a control volume bounded by surfaces of constant buoyancy, we can construct an appropriate background buoyancy profile $b_\ast(z,t)$ and accurately quantify diapycnal mixing in such systems. This new framework is consistent with the definition of a local APE density, useful for identifying mixing mechanisms. The evolution of APE is analysed in various turbulent stratified flow simulations. We show that the mean dissipation rate of buoyancy variance $\chi$ provides a good approximation to the mean diapycnal mixing rate, even in flows with significant variations in local stratification. When quantifying measures of mixing efficiency in transient flows, we find significant variation depending on whether laminar diffusion of a mean flow is included in the kinetic energy dissipation rate. We discuss how best to interpret these results in the context of quantifying diapycnal diffusivity in real oceanographic flows.


Introduction
The transport of heat and salt across surfaces of constant density (isopycnals) in the ocean provides a vital contribution to the closure of the ocean's energy budget (Wunsch & Ferrari 2004;Hughes et al. 2009). As originally highlighted by Munk (1966), such a diapycnal flux arising from molecular diffusion alone is insufficient to balance the generation of dense water in polar regions and close the global circulation. Turbulence therefore plays an important role in enhancing mixing through the stirring of tracer fields (such as temperature or salinity) and the subsequent generation of small-scale gradients. In the ocean interior, turbulence is often associated with breaking internal waves (MacKinnon et al. 2017), which in turn lead to mixing that is strongly intermittent in both space and time. Identifying the mechanisms by which turbulence is generated, and how much mixing can be associated with them, is vital in understanding and accurately modelling the transport of tracers through the ocean.
Here we define mixing as the irreversible diffusive flux across isopycnals that arises due to macroscopic fluid motions, as in Peltier & Caulfield (2003). This flux is sometimes expressed as the mean vertical flux of density, or equivalently buoyancy b = −g(ρ − ρ 0 )/ρ 0 . The flux w b can however include significant contributions from entirely reversible processes such as internal waves. Indeed equating buoyancy flux and irreversible mixing is only appropriate when both are averaged over time and applied to a statistically stationary state. Winters et al. (1995) show that the true rate of irreversible mixing in a Boussinesq fluid is equal to the conversion rate of available potential energy (APE) to background potential energy (BPE). As introduced by Lorenz (1955), APE refers to the change in energy resulting from adiabatically sorting the buoyancy field of a fluid to its state of minimum potential energy. By extending the APE framework to compressible flows Tailleux (2009) argues that mixing should in fact be described as the dissipation of APE into internal energy, which is balanced exactly by an enhancement in the generation of BPE in the Boussinesq limit. In this study, we focus on the dynamics of a single-component Boussinesq fluid with a linear equation of state, and refer the reader to Tailleux (2013a) for a discussion of mixing and APE in more complex scenarios.
Although the Winters et al. (1995) framework provides an exact expression with which to calculate diapycnal mixing, it is not practical for use in oceanographic observations. The most precise observational estimates of mixing come from vertical microstructure profilers that record small-scale gradients of velocity or temperature in isolated vertical profiles. The methods of Osborn & Cox (1972) and Osborn (1980), which are derived from mean balances in the buoyancy variance and turbulent kinetic energy equations respectively, can then be used to estimate an effective diapycnal diffusivity K d . This diffusivity is related to the mixing rate through M = K d N 2 where N is some appropriate measure of the buoyancy frequency. Note that N may not be straightforward to identify if there is significant spatio-temporal variability in the flow (Arthur et al. 2017). Both estimation methods are derived on the assumption that the flow is statistically steady and thus that the mixing is well described by some average of the buoyancy flux. The diffusivity K d obtained from these microstructure measurements can then be checked against estimates of diffusivity from tracer release experiments (Ledwell et al. 2000). Understanding how K d varies throughout the ocean is also vital for improving the accuracy of global circulation models, where diapycnal turbulent fluxes are only captured through a prescribed parameterization of K d , such as that of Klymak & Legg (2010).
Accurately quantifying mixing in computational fluid dynamics requires the use of direct numerical simulations (DNS) that resolve down to the dissipative scales of motion. These simulations can then be used to test the assumptions used to derive the above models (as in Taylor et al. 2019), or to quantify the differences in inferred diffusivity arising from the models (Salehipour & Peltier 2015). The need to resolve the smallest scales of motion restricts the Reynolds numbers Re it is possible to attain through DNS, and so massive computational grids are needed to push Re up towards geophysical values. Since the earliest days of simulating turbulence through DNS, triply periodic domains have been used to reduce computational cost (Orszag & Patterson 1972). The lack of fixed boundaries in this setup means that higher values of Re can be obtained. Thin boundary layers do not need to be resolved and highly efficient pseudospectral methods, exploiting the imposed periodicity, can be implemented. Riley et al. (1981) were the first to include a mean density stratification in such a triply periodic setup by decomposing the buoyancy field into a linear profile N 2 0 z and a periodic perturbation θ. This system has since proved popular for studying the dynamics of high Re stratified turbulence (e.g. Staquet & Godeferd 1998;Riley & de Bruyn Kops 2003;Brethouwer et al. 2007). Investigations of mixing in periodic domains, recent examples of which can be found in Maffioli et al. (2016) and Garanaik & Venayagamoorthy (2019), do not however implement the rigorous Winters et al. (1995) framework for quantifying APE, thus identifying explicitly irreversible mixing. It is common instead to describe mixing in terms of the destruction rate of buoyancy variance χ. This approximation is closely related to the methodology underlying the diffusivity estimate of Osborn & Cox (1972).
As we later explore in §4.3, approximating mixing with χ can result in an over/under-estimate depending on whether the most intense turbulence in the flow preferentially samples regions of locally high/low stratification (and thus is associated with different characteristic local values of the buoyancy frequency). It is therefore important to be able to quantify mixing accurately in the periodic system and identify whether such discrepancies can be significant. Since the buoyancy in the system is only defined through its periodic perturbation θ, ambiguity arises in how to construct the background state of minimum potential energy. In §2 we use a simple example to highlight this issue and then provide an extension to the framework of Winters et al. (1995) that resolves the ambiguity in the case of triply periodic domains. §3 gives a brief overview of the numerical simulations we shall use to test the new framework, including the numerical method used. §4 uses the new framework to analyse the simulations, and compares the exact mixing rates to commonly used estimates. Finally, we conclude and discuss these results in §5, with a particular focus on how our findings impact the estimation and parameterization of mixing in the ocean.

Quantifying mixing in triply-periodic domains
We consider the problem of quantifying irreversible mixing in a system governed by the dimensionless Boussinesq equations subject to an imposed, constant, mean stratification. We decompose the buoyancy field as b = z + θ, where b = z represents the buoyancy profile of the imposed mean stratification. Note that b has been non-dimensionalized by L 0 N 2 0 , where L 0 is a typical length scale and N 0 is the mean dimensional buoyancy frequency, so the mean buoyancy gradient in the dimensionless system is always equal to one.
We apply periodic boundary conditions in all three directions to the flow velocity u, pressure p and buoyancy perturbation θ. The (dimensionless) domain sizes in the x, y, and z directions are L x , L y , and L z , respectively. The dimensionless parameters in the system are the Reynolds number, Prandtl number and bulk Richardson number, given by where U 0 is a velocity scale, ν is the kinematic viscosity, and κ is the diffusivity of buoyancy. As mentioned in the introduction, these equations are frequently used in studies of stratified turbulence where the periodicity allows for the use of efficient spectral methods and removes the effect of solid boundaries. Although the buoyancy perturbation θ is periodic in the vertical, the total buoyancy b = z + θ is not. We are instead left with a jump condition for b at the upper and lower boundaries that has consequences for the calculation of irreversible mixing and potential energy. The mean potential energy in the domain is defined as where f = 1 V ∫ V f dV denotes an average over the domain volume V. Substituting θ = b − z into (2.3) and multiplying by −Ri 0 z provides an evolution equation for the potential energy in the (2.6) Taking the top and bottom boundaries to be at z = L z and z = 0 respectively, and applying the periodic boundary conditions simplifies the equation to where an overbar denotes a horizontal average, defined as f = 1 A ∬ A f dA where A is the crosssectional area of the domain and dA = dxdy is the area element. The conversion rate of internal energy to potential energy, given by the third term on the right hand side of (2.6), has been cancelled out by the main contribution of the diffusive flux through the boundary -the final term in (2.6). The evolution equation (2.7) highlights how sensitive the evolution of the potential energy can be to the choice of the boundary.
The accurate quantification of irreversible mixing requires partitioning the potential energy into background and available components. The background potential energy (BPE) is defined as the minimum value of potential energy that can be achieved through adiabatic rearrangement of the fluid in the domain. In this minimum state, the buoyancy profile is given by a monotonically increasing one-dimensional function b * (z, t), so the mean BPE is given by Winters et al. (1995) show that BPE can also be expressed as where z * is the height a parcel of fluid with buoyancy b(x, t) is moved to under the adiabatic rearrangement. Following Lorenz (1955), the available potential energy (APE) is defined as the surplus potential energy P A (t) = −Ri 0 b (z − z * ) . (2.9) The rate of irreversible mixing associated with fluid motion is then given by the energy transfer rate from APE to BPE, which takes the form where Z * (b, t) is the inverse function associated with the sorted buoyancy profile b * which satisfies z * (x, t) = Z * (b(x, t), t). It is important to appreciate that the term scaling |∇b| 2 in (2.10) is effectively the inverse square of the buoyancy frequency of the sorted variables, and so accentuates the contributions where the sorted buoyancy gradient is relatively weak. As discussed below in §4.3, this is a potential source of difference between M and the buoyancy variance destruction rate χ. Note that Tailleux (2013a) argues for a more general definition of APE, where b * is replaced by an arbitrary reference state b r that can depend on a wide range of thermodynamic quantities. This definition is particularly useful for its possible extension to multicomponent, compressible fluids as shown by Tailleux (2018). However for the single-component, Boussinesq, linear equation of state considered here, we shall continue to use the adiabatically sorted b * given its aforementioned connection to a widely-studied approach used for the quantification of diapycnal mixing.
We now present a simple example to highlight how the aperiodicity of b can cause issues for calculating the mixing rate M. We consider the buoyancy field given by θ = sin x in a domain of length 2π. This might be thought of as a representation of the buoyancy field associated with Note that an overbar here denotes an average over x, and ∂ Z * /∂b is evaluated at b(x, z). a standing internal gravity wave, at an instant when half the columns of fluid in the domain are raised up and half are pushed down relative to their equilibrium location.
The total buoyancy field b = z + sin x and its corresponding sorted profile Z * (b) are shown in figures 1a and 1b respectively. In an unbounded domain, we would expect a linear profile for Z * since the wave is simply a rearrangement of the initial linear stratification. However by taking the boundaries at z = 0 and z = 2π, we produce a profile with deviations from a uniform slope close to these values. Since θ is independent of z, we would also expect the mixing rate M to be constant regardless of the vertical extent that we average over. Figure 1c instead shows that the variations in ∂ Z * /∂b change the value of M across much of the domain, with the horizontally-averaged mixing rate even taking negative values close to the boundary. This issue has caused problems in the literature before. By sorting the buoyancy profile of a rectangular, periodic domain, Bouruet-Aubertot et al. (2001) observe extremely large oscillations in the APE of a breaking internal wave (as shown in figure 9a of that paper) due to fluxes across the top/bottom boundary of the domain. The large oscillations make it difficult to draw precise, quantitative conclusions about the evolution of APE and diapycnal mixing in that system.

Potential energy between isopycnal boundaries
We propose the use of a control volume bounded by surfaces of constant buoyancy (isopycnals) to tackle the highlighted issue of quantifying mixing in triply-periodic domains. Consider tiling the computational domain by stacking several computational domains vertically, as in figure 2. The velocity and buoyancy perturbation repeat in each domain, but the vertical coordinate, z, is continuous such that the total buoyancy in one tile is L z larger than the total buoyancy at the same relative position in the tile immediately below it. In this system it is particularly useful to consider two isopycnals separated by the vertical periodic length, i.e. L z . These isopycnals will have the same shape due to the periodicity of θ, and the volume enclosed by these two isopycnals will therefore be constant. The buoyancy profile can then be sorted into a background state b * (z), where the parcels are sorted into the one-dimensional domain 0 z < L z . Although this background profile must have a mean vertical gradient equal to the imposed mean stratification, its local gradients ∂b * /∂z can vary more generally. In the simple example considered in figure 1, this technique recovers the linear profile Z * (b) = b expected from the column displacement argument mentioned above. , and buoyancy perturbation θ, shown in (c), simply repeat thanks to their periodic boundary conditions. The total buoyancy b = z + θ, shown in (d), is not periodic in the vertical, although isopycnal surfaces separated by the vertical period L z are of identical shape.
We now describe more precisely the details of implementing isopycnal boundaries for quantifying available potential energy and mixing. We first choose a buoyancy value b 0 that defines the lower boundary surface z 1 (x, y, t) implicitly through b(x, y, z 1 (x, y, t), t) = b 0 . (2.11) Vertical periodicity of θ then requires that the upper boundary surface b = b 0 + L z is defined by z 2 = L z + z 1 . It is important to appreciate that (2.11) defines z 1 (and hence also z 2 ) as a single surface that spans the horizontal cross-section of the domain. This ensures that the volume enclosed by the isopycnals is clearly defined. To aid the calculation of volume integrals, we also require (essentially for clarity of exposition) z 1 to be a single-valued function of x and y, or equivalently that the boundary isopycnal cannot exhibit overturning. Such an isopycnal may be difficult to find in homogeneous turbulence, although stratified flows are often strongly spatially inhomogeneous. A discussion of how this approach could be generalised for an overturning isopycnal surface can be found in appendix A. Constructing evolution equations for mean energy quantities involves taking time derivatives of volume integrals. Since the boundaries of our domain are now time-dependent, we must apply the Leibniz rule to any such integral, that is d dt where A is the horizontal cross-sectional area of the domain and the area element dA = dxdy. The mean kinetic energy of the system K = |u| 2 /2 is unaffected by the change of boundaries, since its integrand is periodic in the vertical direction. The evolution of K can therefore be derived straightforwardly from (2.2) by applying (2.1) to obtain the simple form where the buoyancy flux and kinetic energy dissipation rate are respectively given by (2.14) Note that from this definition, positive values of buoyancy flux correspond to a conversion of potential energy to kinetic energy. However, extra terms do arise compared to (2.6) when deriving the potential energy evolution equation. These new terms provide a secondary reservoir of potential energy for the system, as is explained below. The advective flux across the boundary, given by the second term on the right of (2.6), is now zero since the bounding isopycnals have the same shape and the same gradients due to periodicity. Applying the Leibniz result (2.12) to the potential energy P and imposing the boundary conditions therefore produces the evolution equation (2.15) (A more detailed derivation of this equation can be found in appendix B.1.) D p = Ri 0 /RePr is the conversion rate of internal energy to potential energy, and F d is the diffusive boundary term given by where the overbar denotes a cross-sectional average over x and y, importantly taken after the quantity in brackets is evaluated at z = z 1 (x, y, t). We refer to the quantity S as the surface potential energy, where S is defined as (2.17) We can arbitrarily set b 0 = 0 in all of the above by shifting our vertical coordinate to z − b 0 . S then takes the form of potential energy associated with an interface at z 2 , motivating our choice for its name.

APE and BPE between isopycnal boundaries
Using the Winters et al. (1995) form of APE defined in (2.9) is not appropriate for the timevarying domains considered here. This can be understood by considering the simple two-layer system shown in figure 3. Panel (a) shows the background state obtained through constructing the one-dimensional buoyancy profile b * for the buoyancy fields in panels (b) and (c). Since the buoyancy field in figure 3b can be obtained from the background state through shifting the same number of fluid columns up as down, P does not change between states (a) and (b). P = P B therefore holds for state (b), and hence P A = 0. It is simple however to construct a state (c) with lower potential energy than state (b). The Winters et al. (1995) definition would then in fact give P A < 0 for the buoyancy profile in figure 3c, which is not consistent with the concept of available potential energy.
We aim to define a new APE variable A that can be used in the time-varying domain. Progress can be made by considering the total potential energy P + S that appears in (2.15). The decrease in P from figure 3a to 3c is matched exactly by an increase in S. In terms of the total potential energy, states (a) and (c) are therefore equivalent background states. This motivates subdividing the potential energy into P + S = A + B. any change in P + S due to a vertical shift of the domain is captured by the background potential energy B. We therefore construct the background profile b * (z) over the domain z 1 < z < z 2 , such that This ensures that any change in P due to a shift in the mean height of the lower isopycnal z 1 leads to a corresponding change in P B . Accounting also for the corresponding change in S leads to the following definitions for background and available potential energy: Note that for a closed system with fixed, insulated boundaries, these definitions recover the Winters et al. (1995) form for BPE and APE up to a constant in the BPE. Evolution equations for these quantities can be readily obtained through multiplying the buoyancy evolution equation (2.3) by z * and taking volume averages. An analogous derivation as that leading to (2.15), as shown in appendix B.2, results in where M is the irreversible mixing rate defined in (2.10). Subtracting (2.22) from (2.15) also gives an evolution equation for our new APE variable as We therefore recover the simple evolution equation for APE in a closed system, where the irreversible mixing rate M may also be identified with a destruction of APE (e.g. Peltier & Caulfield 2003).

Comparison to local APE of Scotti & White (2014)
The concept of local APE is used as an alternative framework for quantifying available potential energy in situations where fluxes through a boundary are important. Originally devised by Holliday & McIntyre (1981) and Andrews (1981), local APE has seen renewed interest recently in its application to numerical simulations. We follow Scotti & White (2014) in defining the local APE density E APE as a function of space and time by ( 2.24) We use this form primarily for its ease of notation, although as we show in appendix C, for the setup we consider it is equivalent to various other expressions proposed for local APE density. Although this quantity varies in space and time, its dependence on the globally sorted profiles b * and Z * means that it cannot be calculated solely from local fields. In particular, the issue for quantifying mixing highlighted by figure 1 remains unless we change the domain over which b * is constructed. If we instead take isopycnal boundaries as in §2.1, then E APE becomes a useful tool for investigating the local mechanisms that lead to mixing in the domain. Indeed, we can relate the volume-averaged E APE to our global APE variable A as follows. The mean local APE defined in this way can also be written in the form ( 2.25) Winters & Barkan (2013) explain that the final term in this expression accounts for the energy changes arising from the requirement of incompressibility, leading to displacement of some fluid elements to make room for the rearrangement of a fluid parcel in the sorting process. They also showed through considering fluid parcel exchanges that this term vanishes in the case of fixed horizontal boundaries. We now consider a simple example to show how this term can change with non-horizontal boundaries and how it relates to the additional terms in (2.21). We take θ = −z 1 (x, y, t) as the buoyancy perturbation field, so the domain represents that of a uniform stratification where each fluid column has been shifted so that the b = 0 isopycnal is at z 1 , analogously to the situations shown in figures 1a and 3b. In this case, the reference profiles simply take the form b * (s, t) = s − z 1 (t), and Z * (s, t) = s + z 1 (t). We can therefore analytically compute Considering each of the terms in (2.25) separately, we also find that The integral term therefore accounts for all of the available potential energy associated with the surface potential energy S. Since z 2 2 − z 2 2 = z 1 2 − z 1 2 , we find that the volume-integrated local APE exactly matches the changes we propose for the global APE in (2.21). In contrast to A, the calculation of E APE does not rely on computing a surface integral over the isopycnal boundary. Indeed the moving boundary only affects the calculation of local APE through the boundary conditions (2.19) for the reference profiles b * (z, t) and Z * (s, t), and even here one only needs to know the mean height of the isopycnal. The strong agreement between A and E A gives us hope that in flows where A is not well defined, E A can provide an accurate measure of available potential energy.

Numerical simulations
We apply the extended APE framework developed in §2.2 to two sets of direct numerical simulations. All of these simulations are performed using D , which uses a third-order Runge-Kutta scheme for time stepping and a pseudo-spectral method for calculating spatial derivatives (Taylor 2008). The software also implements dealiasing of nonlinear terms through a 2/3 rule. One set of simulations (set F) adds forcing terms to (2.2) and (2.3) to produce a flow in a statistically steady state, whereas the other simulations (in set U) solve the equations unforced Reynolds number (Re) 10 000 10 000 10 000 8000 8000 5000 5000 as an initial value problem. All of the simulations considered here are performed using a bulk Richardson number Ri 0 = 1 and Prandtl number Pr = 1. The first set of simulations are those used in our previous study on mixing in forced stratified turbulence (Howland et al. 2020). We refer to the simulations H, R and P in that paper by F1, F2, and F3 respectively, and outline some of their key parameters in table 3. Simulation F1 is forced by randomly phased large-scale vortical modes, and importantly features no direct forcing of the buoyancy field. The evolution equations (2.22) and (2.23) for B and A still therefore hold. On the other hand, simulations F2 and F3 are forced by large-scale internal gravity waves that include a buoyancy forcing component. The buoyancy forcing can act as a source or sink of potential energy, modifying the evolution equations. However if we are primarily concerned with diapycnal mixing, it remains useful to calculate the irreversible mixing rate M in these cases. For more precise details of the forcing in these simulations, we refer the reader to Howland et al. (2020).
The second set of simulations investigate the interaction of a sinusoidal vertical shear flow and a plane internal gravity wave. The initial velocity and buoyancy fields are given by u = (sin z) x + u and θ = θ respectively, where We express the initial amplitude of the internal wave through its steepness s and choose the wave vector k = (k, l, m) = (1/4, 0, 3) based on the typical aspect ratios of waves observed in the thermocline by Alford & Pinkel (2000). Small-amplitude noise is added to the initial velocity field to allow the development of three-dimensional motion from the two-dimensional initial condition. Simulations U1 and U3 use an initial wave steepness of s = 1, with s = 0.5 for simulation U2 and s = 0.75 for simulation U4. As an example, the initial buoyancy and spanwise vorticity fields for U4 are shown in

Energy budgets
We now investigate the evolution of background and available potential energy in the various simulations, and consider how terms in the energy budgets (2.22) and (2.23) relate to the flow dynamics. Figure 5  (2.22) are close to being equal, as shown in both figures 5 and 6. The diffusive boundary term F d primarily acts to cancel out any increase in B due to the conversion of internal energy to potential energy through D p . This cancellation is exact when the boundary has no lateral variation, and arises since the system is forced to maintain a constant mean buoyancy gradient through the periodicity of θ. Figure 6 repeats the analysis of figure 5 for the forced simulation F1. We only consider the statistically steady period achieved at late times in this flow. Unlike in the unforced simulation, the mean buoyancy flux remains negative throughout as shown in figure 6a, providing a source of APE from the kinetic energy. Figure 6b furthermore shows that the buoyancy flux is on average in balance with the mixing rate, leading to an approximately constant value of A, as shown in figure 6c. The constant mixing rate also predictably leads to a linear increase in the background potential energy.

Visualising mixing with local APE
We can further investigate the local processes that lead to the global results above by analysing the distribution of local APE throughout the domain. Figure 7 plots snapshots of E APE (x, t) from simulations F1-F3 and from simulation U1 at various times. Since the turbulence arising in each simulation is patchy and inhomogeneous, we are able to choose appropriate isopycnal boundaries for each simulation and hence calculate the surface potential energy S. These isopycnal boundaries are shown in figure 7 as solid black lines.
Data from the forced simulations of set F are presented in figures 7a-c. Each snapshot of E APE is taken at time t ≈ 150, when the turbulence is in a statistically steady state. Figure 7a highlights low local APE values throughout the domain of simulation F1. Increased E APE occurs only at small scales and in regions with high turbulent dissipation rates (not shown). In this sense, APE is primarily associated here with the distortion of the buoyancy field by turbulence, and not with internal waves. By contrast, figures 7b and 7c show patches of high local APE throughout the domain at a range of scales. This is consistent with associating mixing with intermittent, large-scale overturns and convectively-driven turbulence, as discussed in Howland et al. (2020).
The development of local APE during the unforced simulation U1 is presented in figures 7d-f. The distribution of E APE in the initial condition is shown in figure 7d, and is entirely associated with the internal gravity wave described by (3.1). At early times, the wave is refracted by the shear flow, leading to a distortion of the banded structure in the local APE field. By time t = 20, E APE preferentially accumulates in the upper half of the domain while maintaining some signal of the wave structure, as shown in figure 7e. The large values of E APE lead to locally unstable buoyancy profiles, and the development of convective instabilities. The associated convection converts APE to kinetic energy through the buoyancy flux, and also promotes the emergence of small scale structures seen in figure 7e. Later, at t = 30, the flow becomes more complex with the development of shear-driven turbulent billow structures. These structures, seen prominently on the right of figure 7f, span regions of both high and low E APE . Although the volume-averaged mixing rate peaks near this time, the banded structure of E APE leads to strong local variation in local mixing rates within the turbulent patches. Mixing is high where turbulence and APE coexist, and it cannot occur where there is no APE to remove.

Estimating mixing with χ
In the limit of small buoyancy perturbations from the uniform, imposed buoyancy gradient, available potential energy can be approximated bỹ This quantity satisfies the simple evolution equation where χ is the rate of destruction of buoyancy variance, i.e. the dissipation rate defined by Comparing the definitions (2.10) and (4.3), we see that χ precisely takes the form of the irreversible mixing rate M only for the specific case where the sorted background buoyancy profile b * exactly matches the imposed uniform stratification. Recall that this imposed constant stratification has a dimensionless buoyancy gradient equal to one by construction. In our simulations, local deviations in the buoyancy field are not always small and we should treat the above approximation with caution. For example, during the convective phase (20 < t < 30) of simulation U1 there are sizeable regions of the domain with statically unstable buoyancy gradients. The peak mixing in this case occurs where the (horizontal) mean buoyancy is in a layered state, with 'layers' of relatively low stratification separated by 'interfaces' of relatively high stratification compared to the imposed constant buoyancy gradient. Such layered states are observed to arise naturally in turbulent stratified flows, for a wide variety of dynamical reasons (see Caulfield 2021, for a review). Nevertheless the dissipation rate χ is significantly more straightforward to quantify than the true mixing rate M, and so it is useful to investigate how well it can actually approximate the mixing. The accuracy of χ for estimating mixing is also important in the context of ocean microstructure measurements, where small-scale gradients are measured directly but there is no way to obtain the relevant reference profile b * . In figures 8a-c, we therefore plot the time series of both χ and M for each of our simulations. By inspection, the two quantities appear to match up very well, with the symbols marking the mixing rate overlapping the lines plotting the time series of χ. To quantify how well χ approximates M, we plot the time series of their ratio in figures 8d-f. Throughout the forced simulations, and for the early times of the unforced simulations, χ remains within 10% of the true mixing rate. At late times in simulations U1 and U3 the difference increases up to 20%, but at this stage the flow is relaminarizing and M and χ are both small. Indeed we show that this discrepancy is unimportant for quantifying the total mixing achieved over the course of the simulations in figures 8g-i, where we plot the ratio of time-integrated χ and M. The time integral of M is equal to the increase in background potential energy due to diapycnal mixing, and we see that using χ to estimate this quantity results in at most a 5% error in the total BPE change (corresponding to the final values of the cumulative ratio plotted in figures 8g-i).
In the unforced simulations, χ consistently provides a slight underestimate of the diapycnal mixing rate. This suggests that regions of intense turbulent mixing, associated with high values of |∇b|, preferentially sample regions where ∂ Z * /∂b > 1. These regions are in turn associated with the reference buoyancy profile b * having a locally weaker stratification than the mean. In simulations F2 and F3, where forcing is applied in the form of internal gravity waves, the opposite is true and χ provides a slight overestimate of M. However it is not true that intense mixing occurs only in regions of strong or weak local stratification in each flow. In all of the forced simulations for example, the standard deviation of ∂ Z * /∂b rises from the range 0.1-0.15 at time t = 50 up to 0.25-0.3 at t = 150, suggesting that as mixing persists throughout the simulations, the background profile is modified. The fractional error between χ and M seen in figure 8d does not show this increasing trend, suggesting that some local overestimates of M (where ∂ Z * /∂b < 1) cancel with some local underestimates (where ∂ Z * /∂b > 1) in the global average. Similarly, the standard deviation of ∂ Z * /∂b reaches values in the range 0.15-0.2 for simulations U1 and U3 when t > 30, approximately double the fractional error during the period of peak mixing.

The effect of mean flow dissipation
In the unforced simulations of set U, the majority of the kinetic energy is associated with the initial mean shear profile u = sin z. At late times in these scenarios, the flow begins to relaminarize and the kinetic energy dissipation rate ε is dominated by the laminar diffusion of the mean shear. Mixing efficiency is however often calculated using the turbulent kinetic energy dissipation rate, that we quantify here as where u = u − u is the velocity field perturbation from the horizontal average. Figure 9 compares time series of the following definitions of mixing efficiency calculated using either the turbulent dissipation rate ε or the total rate ε We use χ rather than M in our definition of efficiency, since we have seen that the difference between them is small in the previous section, and our records of χ have better resolution in time.
Large discrepancies between η and η are observed when the average TKE dissipation rate ε is small compared to the dissipation rate of the mean flow ε = |∂ u/∂z| 2 /Re. In simulation U2, wave breaking occurs at t ≈ 50 and consists of small, strongly localised overturns that dissipate relatively quickly. Consequentially ε remains smaller than ε for the entire duration, leading to large differences between the efficiencies in figure 9b. η takes much larger values than η in all of the unforced simulations at early and late times, with η close to its initial value of 0.5. This value corresponds to the diffusion associated with the plane wave form of (3.1) and is a consequence of the choice Pr = 1, that is molecular diffusion of buoyancy occurs at the same rate as the diffusion of momentum. At larger values of Pr, diffusion of the wave would result in a far lower value of η .
In figure 9d-f we also plot associated cumulative mixing efficiencies, defined here in terms of appropriate integrals of χ and ε (or ε ): where t 0 = 50 for the forced cases, and t 0 = 0 for the unforced cases. The time integrals represent the energy changes associated with the cumulative effects of χ and ε. Figures 9e and 9f show that the diffusion of the mean shear flow has a significant impact on the total cumulative efficiency in the unforced simulations. To emphasise that this is primarily a Reynolds number effect, we have listed measures of the buoyancy Reynolds number for the unforced simulations in table 4.4.
Since the flows considered are extremely inhomogeneous in the vertical (as seen in figures 7e and 7f) we have calculated Re b from both volume and horizontal averages. In oceanographic flows, we expect molecular diffusion to be negligible compared to the turbulent dissipation rate ε Re Ri 0 11.29 0.46 6.89 2.41 for the vast majority of the internal wave spectrum. This result therefore highlights the challenge of using direct numerical simulations, where Re is inevitably limited by computational resources, to investigate ocean mixing processes.

Discussion and conclusions
In this study, we have highlighted how the APE framework of Winters et al. (1995) should be generalised in the triply-periodic system often used in numerical simulations of stratified turbulent flows. In these systems it is important to constrain the buoyancy field, inferred from the periodic perturbation θ, to lie in a prescribed range. We can then construct an accurate background buoyancy profile b * that is consistent with the periodic nature of the system. However, setting limits on the buoyancy values effectively means that the shape of the domain can change in time. In the case where the limiting buoyancy value has a non-overturning isopycnal surface, we find that this introduces an extra potential energy term S as defined in (2.17). Appropriate definitions of available and background potential energy can then be obtained by accounting for this additional term as in (2.20) and (2.21).
Constructing the correct background profile is also vital for accurately calculating the local APE density E APE defined by Scotti & White (2014). This quantity can then provide useful information for identifying mechanisms by which mixing can occur. When integrated over the domain, the local APE also recovers all of the additional terms in our new global APE variable A. Furthermore, the local APE can even be quantified in scenarios where our global APE is not well defined. So long as the background profile b * is identified, both E A ≡ E APE and the irreversible mixing rate M can be calculated. The evolution of E A is then entirely determined by the mixing rate and the buoyancy flux, with zero contribution from boundary fluxes. We can therefore calculate the exact rate of diapycnal mixing in more energetic stratified flows that use periodic domains, such as those considered by de Bruyn Kops & Riley (2019) and Portwood et al. (2019).
This technique for calculating APE could also be applied to unstably stratified periodic systems, where Ri 0 < 0, used to study bulk properties of convection (e.g. Lohse & Toschi 2003). In traditional Rayleigh-Bénard convection, Gayen et al. (2013) find that irreversible mixing is largely confined to thermal boundary layers. It would therefore be interesting to investigate whether the theoretical prediction of η → 0.5 at high Ra holds for the periodic convection setup, where such boundary layers are absent. The sorting technique presented here is also applicable to the case of passive scalar flows where a mean gradient is imposed. Such setups are useful for studying the mixing of biogeochemical tracers, which are often found with significant mean gradients in the ocean (Williams & Follows 2011). Although the concept of APE would not be relevant here, sorting between isoscalar surfaces would provide an appropriate background profile to enable accurate calculation of the diascalar flux as proposed by Winters & D'Asaro (1996).
In observational oceanography, turbulent mixing can be estimated by using fast-response thermistors to measure small-scale temperature gradients. The primary aim in this context is to estimate a diapycnal diffusivity, defined in our dimensionless formulation as Since the reference profile b * cannot be obtained in the ocean, a large-scale average is taken of the buoyancy (or temperature) gradient. The estimate often attributed to Osborn & Cox (1972) is then used such that Note that the internal energy conversion rate D p is neglected here, since it is assumed to be much smaller than χ in a turbulent flow. In dimensional form it is common to see (5.2) written as K d = χ/N 2 , but in our non-dimensionalization the mean buoyancy gradient in the denominator is prescribed to be equal to one. The approximation made in estimating K d in (5.2) is the same approximation used in §4.3 to estimate the mixing rate M with χ. Precisely, we approximate the reference buoyancy gradient ∂b * /∂z by the imposed mean stratification. We test this approximation in the context of diapycnal diffusivity in figure 10 by plotting the time series of ( χ + D p )/K d . The fractional error between the estimate χ + D p and the true diffusivity remains within one standard deviation of ∂ Z * /∂b for every simulation. Figures 10b  and 10c show that χ + D p underestimates the diffusivity at the time of most intense mixing in the unforced simulations. This reaffirms the conclusion drawn from figures 8e and 8f that the turbulent mixing in this flow preferentially samples regions of relatively weak local stratification.
Salehipour & Peltier (2015) find a similar underestimation of K d in turbulent flows developing from Kelvin-Helmholtz instability in a stratified shear layer. An investigation to identify in which flows (5.2) provides an over/underestimate of the diffusivity would be valuable for understanding the variability associated with the approximation. We include the internal energy conversion rate D p in our estimate in figure 10 since it is not always negligible in the simulations. Furthermore Gregg et al. (2018) remark that D p should be included when applying mixing results to the strongly stratified pycnocline where mixing is localised and intermittent. In the periodic setup studied here, the boundary flux F d counteracts D p in the BPE energy budget (2.22) to maintain the constant mean stratification. When quantifying diffusivity in this system, it is therefore important to include the contribution from D p and to compute M + D p directly, instead of relying on changes in BPE. In many observational studies focused on mixing in turbulent patches (where D p is negligible), practical difficulties in obtaining an accurate value of χ result in far larger implied levels of uncertainty than those apparent in in figure 10 (see for example Waterhouse et al. 2014). In this sense our results show that (5.2) provides a good estimate of the diapycnal diffusivity in the stratified flows considered.
In the case of homogeneous turbulence subject to a uniform mean stratification, Stretch & Venayagamoorthy (2010) show χ and M to be equivalent. Indeed, if such homogeneous turbulence is maintained in a steady state by energy transfers from the velocity field, then χ is also equivalent to −J . This is exactly the reasoning of Osborn & Cox (1972). In the periodic system considered here, if mixing were homogeneous throughout the domain then the boundary flux F d would also balance the interior diapycnal flux M + D p such that the BPE equation (2.22) was steady. As highlighted by Portwood et al. (2016), many stratified flows are not homogeneous in this sense, with turbulence becoming more patchy and intermittent when subject to stronger stratification. This even applies to flows where the initial state is homogeneous, such as the decay of a turbulent cloud in an initially uniform stratification (e.g. Bartello & Tobias 2013). We believe the APE framework presented above will prove useful in determining the potential impacts of such developing inhomogeneities.
Due to the aforementioned difficulties involved in accurately resolving small-scale temperature gradients, shear probes are used more frequently than thermistors to infer mixing rates in the ocean. Further assumptions are however needed to obtain mixing estimates from such velocity gradient measurements. On top of the Osborn & Cox (1972) model, the buoyancy variance destruction rate may be approximated by χ −J = Γε, where the turbulent flux coefficient Γ is taken to be a constant, usually 0.2 in practice after Osborn (1980), under a set of assumptions that the turbulent flow is, for example quasi-steady. The turbulent flux coefficient is related to the mixing efficiency defined in (4.5) through Many experimental and numerical studies have shown variation in the mixing efficiency across a range of stratified flows, as reviewed by Ivey et al. (2008) andCaulfield (2021). This has motivated a body of work to investigate the functional dependence of η on various dimensionless parameters, including the Richardson number, buoyancy Reynolds number Re b = ε/νN 2 , and turbulent Froude number Fr = ε/NK. Despite this concerted effort to provide insight into how η varies, there is no clear physical explanation as to why Γ = 0.2 is a sensible assumption or why it appears to provide diffusivity estimates in line with those from tracer release experiments (Ledwell et al. 2000). In figure 9 we highlight examples where laminar diffusion of a shear flow can strongly impact the calculated values of η. Although not relevant for high Reynolds number flows found in the ocean, it is important to acknowledge the effect of this diffusion in idealised numerical studies that discuss mixing efficiency in the context of ocean mixing. This is most relevant for flows where turbulence is transient and localised, such as those arising from instabilities in stratified shear layers.
In the context of estimating mixing from oceanic measurements, the poorly constrained variability of Γ implies that the Osborn & Cox (1972) model will instead provide the best estimates of mixing from microstructure data. Indeed our results above show that in turbulent flows with a large-scale mean buoyancy gradient N 2 0 , the Osborn & Cox (1972) model (5.2) provides a reliable estimate of the diapycnal diffusivity. In oceanic flows calculating N 0 can however prove challenging, particularly in the case of internal waves breaking close to boundaries as highlighted by Arthur et al. (2017). Given the significant dissipation of energy close to boundaries in the ocean, the calculation of N 0 in such flows remains an important outstanding issue for estimating diapycnal mixing in these regions. Furthermore, the Osborn & Cox (1972) method can only be applied in regions where the ocean is stratified by temperature. Where salt acts as a stratifying agent, the turbulent flux coefficient Γ must be specified to obtain microstructure mixing estimates through the Osborn (1980) method.
In particular, for the energetic framework presented here to be truly applicable to real oceanographic flows, there are at least three open issues which need to be addressed. First, it is not at all clear what the effect of more realistic Reynolds numbers, or indeed realistically higher values of Pr = O(10 − 1000) will have on the various mixing properties and energetic pathways discussed here. Second, it is still an open question of some controversy whether Γ ≈ 0.2, or equivalently η ≈ 1/6, is actually 'typical' of quasi-steady mixing processes, or whether Γ actually depends on parameters of the flow. Portwood et al. (2019) recently demonstrated the emergence of Γ = 0.2 in sheared DNS that was controlled by construction to be quasi-steady. It is at least plausible that the higher values of efficiency observed for the flows discussed here are artefacts of the inherent transience of these flows. Of course, mixing events in the ocean are likely to be highly spatio-temporally intermittent, not least because of the key role played by 'breaking' internal waves, as argued by MacKinnon et al. (2017) and modelled here, so the relevance of quasi-steady sustained stratified turbulence to the real ocean is not immediately obvious. Thirdly, complications associated with layered states, either due to hydrodynamic mechanisms associated with turbulence (Caulfield 2021) or associated with double-diffusive convection (Schmitt 1994) are clearly of interest. The energetic framework presented here is nevertheless well-suited to address these three open issues, or indeed other challenges of real relevance to the quantification and parameterization of mixing in realistic stratified flows.

Declaration of interests
The authors report no conflict of interest.

Appendix A. Consideration of a more general boundary isopycnal
In (2.11), we assume that the boundary buoyancy contour can be parameterized by x and y. Now let us consider a more general isopycnal boundary that may overturn, where the surface of constant buoyancy is parameterized by arbitrary coordinates p and q. The implicit definition of the isopycnal surface x 1 (p, q, t) is then given by b(x 1 (p, q, t), y 1 (p, q, t), z 1 (p, q, t), t) = b 0 .
(A 1) Considering the same volume integral as in (2.12), we apply the Reynolds Transport Theorem to obtain d dt S denotes the domain in (p, q) space that parameterizes the surface, x 1 = (x 1 , y 1 , z 1 ) is the location of the isopycnal surface in Cartesian coordinates, and the area element is given by Note that for p = x and q = y, this recovers the original Leibniz rule result of (2.12) since x 1 = (x, y, z 1 (x, y, t)) and We know in general that the direction of the normal is that of the buoyancy gradient ∇b, but for the arbitrary form (A 3) the magnitude of x p × x q depends on the coordinates chosen. Since we wish to calculate the surface integral from simulation data, it is convenient to restrict ourselves to non-overturning isopycnals, where the magnitude of the area element can be straightforwardly obtained.
We can however manipulate (A 2) further by noting that n = ∇b/|∇b|, and defining the average over the surface S as where A S is the surface area of the isopycnal defined in (A 1). Applying this to the Reynolds transport theorem result (A 2) gives d dt Substituting f = −Ri 0 bz to find the extra term in the potential energy equation provides where A is the cross-sectional area of the domain in the x-y plane, and from Winters & D'Asaro (1996) we know that Although the last term in (A 7) can be expressed analytically, its computation is far more arduous than −dS/dt, and it does not appear (thus far) to simplify to a similar form.

B.1. Total potential energy
In this section, the control volume is bounded by the isopycnals b = b 0 (z 1 ) and b = b 0 + L z (z 2 ). Consider the time evolution of P = −Ri 0 bz by applying the Leibniz rule as in (2.12): We can show that this final integral is zero by considering the evolution of the volume-averaged buoyancy. Since b = z + θ, we know that b = L z /2 + z 1 + θ . The mean buoyancy perturbation is coupled to the mean vertical velocity through the system Importantly, if both θ and w are initially zero, then they remain so forever. This is the scenario most commonly applied in studies using the periodic stratified setup, so we proceed taking θ ≡ 0. We therefore know that Applying the Leibniz rule of (2.12) to b instead gives We can then deduce that the desired integral is zero as follows 0 = ∂b ∂t = −∇ · (bu) + 1 RePr ∇ · ∇b , (B 17) where we have applied the divergence theorem and used that ∇b| z 1 = ∇b| z 2 .

B.2. Background potential energy
In this section, we set b 0 = 0 so the boundary surfaces z 1 and z 2 correspond to the isopycnals b = 0 and b = L z . We begin by determining the time evolution of P B = −Ri 0 bz * . Applying the Leibniz result of (2.12) to this quantity gives The second term in the line above is zero in the case of fixed, insulating, horizontal boundaries. We therefore consider the simple case of θ = −z 1 (x, y, t) to investigate the contribution this term has in the case of time-dependent isopycnal boundaries. As in §2.3, this example has the linear background profiles Z * (s, t) = s + z 1 , and b * (s, t) = s − z 1 , so z * (x, t) = Z * (b(x, t), t) = b(x, t) + z 1 (t) = z − z 1 (x, y, t) + z 1 (t).

(B 27)
Here we have introduced the Casimir that satisfies ∇ψ = z * ∇b. Since Z * is the inverse of b * , and we know that b * (z) = b(x) , we can furthermore deduce that