A general criterion for the release of background potential energy through double diffusion

Double diffusion occurs when the fluid density depends on two components that diffuse at different rates (e.g. heat and salt in the ocean). Double diffusion can lead to an up-gradient buoyancy flux and drive motion at the expense of potential energy. Here, we follow the work of Lorenz (Tellus, vol. 7 (no. 2), 1955, pp. 157–167) and Winters et al. (J. Fluid Mech., vol. 289, 1995, pp. 115–128) for a single-component fluid and define the background potential energy (BPE) as the energy associated with an adiabatically sorted density field and derive its budget for a double-diffusive fluid. We find that double diffusion can convert BPE into available potential energy (APE), unlike in a single-component fluid, where the transfer of APE to BPE is irreversible. We also derive an evolution equation for the sorted buoyancy in a double-diffusive fluid, extending the work of Winters & D’Asaro (J. Fluid Mech., vol. 317, 1996, pp. 179–193) and Nakamura (J. Atmos. Sci., vol. 53 (no. 11), 1996, pp. 1524–1537). The criterion we develop for a release of BPE can be used to analyse the energetics of mixing and double diffusion in the ocean and other multiple-component fluids, and we illustrate its application using two-dimensional simulations of salt fingering.


Introduction
Double diffusion occurs when the density of a fluid is a function of two scalars that diffuse at different rates. Our primary motivation is double diffusion in the ocean, so we will refer to the faster-diffusing scalar as temperature and to the slower-diffusing scalar as salt. Double diffusion can drive a variety of flows, including diffusive convection, salt fingering, thermohaline intrusions and thermohaline staircases (Turner 1985;Schmitt 1994;Radko 2013), and may impact canonical flows, including gravity currents, plumes and Kelvin-Helmholtz billows (Smyth, Nash & Moum 2005;Konopliv & Meiburg 2016;Dadonau, Partridge & Linden 2020). Here, we will analyse the energetics of double-diffusive fluids using the concepts of 'background' and 'available' potential energy.
The background potential energy (BPE) can be found by adiabatically sorting the density field into a monotonically decreasing function of height (Lorenz 1955). The BPE is then defined as the potential energy associated with the sorted density field, and the available potential energy (APE) is the difference between the potential energy of the unsorted density and the BPE. The budgets for APE and BPE were first derived by Winters et al. (1995) for a single-component, incompressible, Boussinesq fluid. Winters et al. (1995) showed that energy can be transferred from APE to BPE, but not vice versa (and hence termed irreversible mixing), via a diapycnal (across surfaces of constant density) diffusive buoyancy flux. This framework was used to formalize the definition of mixing efficiency, an often-used concept in ocean mixing (e.g. Peltier & Caulfield 2003;Gregg et al. 2018). Here, we extend the framework from Winters et al. (1995) to include double-diffusive effects.
A local definition of APE was introduced by Tailleux (2013) and Scotti & White (2014), which was further generalized to include compressibility, a nonlinear equation of state and an arbitrary number of scalar components (Tailleux 2018b). Recently, Tailleux (2018a) derived a new expression for the local APE dissipation rate in a binary compressible fluid with a nonlinear equation of state and showed that in this case the APE dissipation is irreversible. Tailleux (2013Tailleux ( , 2018a included terms to represent diabatic heat and salt fluxes, but did not explicitly consider diffusion or double-diffusive effects. Here we will follow the derivation in Winters et al. (1995) and extend their analysis of the global APE and BPE budgets to include double diffusion. Smyth et al. (2005) considered the problem of double diffusion by applying the Winters et al. (1995) formulation to temperature and salinity separately, defining a background potential energy for each component by sorting temperature and salinity independently. While this approach is useful for quantifying irreversible mixing for each scalar, temperature and salinity do not have a gravitational potential energy separate from the fluid density. This approach also does not capture the single-component limit of equal molecular diffusivities. For example, consider a situation where non-zero horizontal temperature and salinity gradients are compensating such that the density depends only on the vertical coordinate. If the molecular diffusivities of each component are equal and the equation of state is linear, then the analysis of Winters et al. can be applied to the fluid density; if the density profile is stable, the APE will be zero. However, the APE calculated from the temperature and salinity individually will be non-zero.
Here we offer three main results. First, we consider the general criterion for the diapycnal buoyancy flux in a double-diffusive fluid and show that its sign can be expressed in terms of what we call the 'gradient ratio', G ρ ≡ (α|∇T|)/(β|∇S|), as well as the angle θ made between the scalar gradients (figure 1), and the ratio of the molecular diffusivities (κ T , κ S ). Second, we extend the work of Winters et al. (1995) by deriving the volume-averaged APE and BPE budgets for an incompressible Boussinesq stratified fluid with a linear equation of state and two diffusing components. We show that the criterion for an up-gradient buoyancy flux obtained earlier can be related to a condition for a transfer of BPE to APE. Note that, while this result changes the interpretation of the partition between APE and BPE, we will still use the terms APE and BPE to refer to the standard definitions based  on the sorted buoyancy. The criterion that we derive is useful for identifying when double diffusion qualitatively changes the energy transfers in a multiple-component fluid. Third, we generalize the evolution equation for the sorted buoyancy profile (Nakamura 1996;Winters & D'Asaro 1996) for a double-diffusive fluid and relate an up-gradient diapycnal buoyancy flux to a negative effective diffusivity in sorted buoyancy coordinates. We then test our criteria using a simulation of salt fingering (figure 2) and show that the previously observed negative diffusivity of salt fingering (e.g. St. Laurent & Schmitt 1999) can be described in terms of G ρ and θ, generalizing the previous results of Veronis (1965) and Garrett (1982).

Governing equations
We will consider the incompressible Boussinesq Navier-Stokes equations where D/DT = ∂/∂t + u · ∇ denotes the material derivative and u = (u, v, w) is the velocity with respect to Cartesian coordinates (x, y, z). The fluid pressure is p,k is the unit vector in the z direction, b is the buoyancy and ∇ · τ is the divergence of the viscous stress tensor. Additionally, we consider two advection-diffusion equations for the two stratifying elements, and a linear equation of state that relates these quantities to the fluid density, ρ, and hence buoyancy, b ≡ −g(ρ − ρ 0 )/ρ 0 , where g is the gravitational acceleration and ρ 0 is a constant reference density. Specifically, these equations are FIGURE 2. Buoyancy (a) and diapycnal buoyancy flux (b,c) from a two-dimensional simulation of salt fingering. Dashed and solid contours in panels (a,b) show temperature and salinity, respectively. Panels (c) are scatterplots of the diapycnal buoyancy flux (coloured as in (b)) in (G ρ , θ) space. The panels to the right of (a,b) show the sorted height (z * ) averages of the panels and the initial profiles are indicated with a dashed line.
Note that T and S can be viewed as generic diffusing scalars, with molecular diffusivities (κ T , κ S ).

2.2.
Condition for an up-gradient buoyancy flux By introducing a new variable b p , the buoyancy evolution equation can be written as (2.6) We will refer to b p as the 'buoyancy flux potential', since ∇b p is the diffusive buoyancy flux. As illustrated in figure 1, ∇b p is not necessarily in the same direction as ∇b (unlike in a single-component fluid), which can result in a buoyancy flux along isopycnals (surfaces of constant density). We can divide the buoyancy flux into diapycnal and isopycnal components, wheren = ∇b/|∇b| is the unit normal to the isopycnal surface, andt is a unit vector tangent to the isopycnal surface, illustrated in figure 1. Introducing the gradient ratio G ρ ≡ (α|∇T|)/(β|∇S|), the diapycnal buoyancy flux can be written as where cos θ = −∇T · ∇S/(|∇T||∇S|) and θ(x, t) is the angle between the vectors ∇T and −∇S. The scalar angle is defined such that, when θ = 0, 2π, the gradients in T and S point in opposing directions such that they contribute to the buoyancy gradient constructively. This definition implies that θ is also the angle made between the T and S isoscalar surfaces (see figure 1). We will refer to the term in brackets in equation ( The function f sets the sign of the diapycnal buoyancy flux, ∇b p ·n, since the other terms in equation (2.8) are positive definite. Therefore, a condition for an up-gradient diapycnal buoyancy flux is f < 0. The dividing line between positive and negative diapycnal flux, f = 0, is a quadratic equation in G ρ and cos θ, plotted in bold in figure 2. The contours of f are also plotted in figure 2. Since G ρ , κ T and κ S are positive, cos θ is the only term in equation (2.9) that can be negative. Therefore, the relative orientation of the surfaces of constant T and S is of central importance to the sign of the diapycnal buoyancy flux. A Reynolds-averaged version of f was derived as a term in the equation for the density variance by de Szoeke (1998), who showed that negative values of this quantity can lead to the production of density variance. The significance of the condition in equation (2.9) to the APE and BPE budgets will be discussed in the next section. When θ = π, the condition for a negative diapycnal buoyancy flux reduces to is the density ratio. This is a well-known result in the double-diffusive literature (e.g. Veronis 1965;Garrett 1982;St. Laurent & Schmitt 1999;Radko 2013). Therefore, equation (2.9) can be viewed as a generalization of equation (2.10), which includes horizontal T and S gradients. Note that salt fingering instability occurs when R ρ > 1 in the case of vertically aligned gradients such that the diffusive buoyancy flux will be down-gradient (Radko 2013). However, as we will demonstrate below, equation (2.9) can be used to analyse the energetics of active salt fingering convection by considering the T and S gradients associated with individual salt fingers. The third row of figure 2 shows contours of f (G ρ , θ ), where the region of f < 0 is represented with dashed contours and the bounding line for this region is in bold. Superimposed are data from a simulation of salt fingers, which will be discussed below. The ∇b p ·n < 0 region is bounded by a critical angle θ c , where (2.11) For κ T = 10 −7 m 2 s −1 and κ S = 10 −9 m 2 s −1 , typical of the ocean, the critical angle is θ c 100 • . In the limit as κ T /κ S → 1, the critical angle θ c → π.
2.3. Potential energy budget In this section we extend the framework of Winters et al. (1995) to include double-diffusive effects. Within a static volume V with boundary S, we can define the kinetic energy (2.12), potential energy (2.13), background potential energy (2.14) and available potential energy (2.15): 14) where z * is the sorted buoyancy coordinate, with H the Heaviside function and A the cross-sectional area of the volume V. The sorted height z * (x, t) can be interpreted as the height of a fluid parcel after sorting the buoyancy field b into a stable configuration, i.e. b(z * ) is the sorted buoyancy profile. Alternatively, z * (b) is the normalized volume of fluid with buoyancy less than b. In practice, z * can be calculated by integrating the probability density function of b (Tseng & Ferziger 2001). Following the method outlined in Winters et al. (1995), the budgets of the energy components for a double-diffusive fluid are (2.20) These budgets can be derived by taking time derivatives of equations (2.12)-(2.15) and applying integration by parts. Additionally, one must use the relations ∇z * = (dz * /db)∇b and dz * /dt z * = 0. By defining ψ ≡ b z * (s) ds, we also use the relation ∇ψ = z * ∇b. These equations are presented schematically in figure 1. Although we follow the derivation in Winters et al. (1995), these equations may also be derived by volume integrating the local APE budgets derived in Tailleux (2018b). The letter S in equations (2.17)-(2.20) denotes a surface flux, with a superscript adv denoting an advective flux, or diff denoting a diffusive flux. The term P adv b is the volume integral of the vertical advective buoyancy flux bw and is the volume-integrated dissipation rate of kinetic energy. The term Φ i is the exchange between internal energy and potential energy. This term is discussed extensively in Konopliv & Meiburg (2016) in the context of the potential energy associated with the horizontally averaged buoyancy. In our volume-averaged formulation, Φ i only involves boundary quantities and in a turbulent flow we might expect Φ i to be small relative to Φ d (Peltier & Caulfield 2003).
Equations (2.17)-(2.20) differ from the single-component case presented in Winters et al. (1995) in the expressions for the diffusive surface terms S diff PE , S diff BPE , S diff APE and Φ i and the APE/BPE exchange term Φ d , where b p now replaces κb. As a result, Φ i can take negative values in a fluid where the density increases with depth, which is not possible in a single-component fluid.
We define the quantity φ d as the integrand of Φ d , i.e.
and using equation (2.8) the sign of φ d is set by the sign of f (G ρ , θ , κ T /κ S ) since dz * /db > 0 by construction. There is a close connection between φ d and the local diapycnal buoyancy flux ∇b p ·n. Start by considering the average of an arbitrary continuous function, g(x, t), over a surface of constant z * , where S * is a surface with constant z * and A s is its area. Next, consider two isopycnals with buoyancy b and b + b and let n and z * denote the perpendicular and sorted distances between the isopycnals. Then, taking the limit as b → 0, we can write where V * = A z * is the volume between the isopycnal surfaces. Winters & D'Asaro (1996) derived this relation for g(x, t) = κ|∇b| (note that their result included a minus sign since they used density instead of buoyancy as the sorted variable).
Taking g(x, t) = ∇b p ·n and using equation (2.23) gives and hence φ d z * is related to the mean diapycnal buoyancy flux. Note that φ d z * and ∇b p · n z * have the same sign since A and A S are positive. Similarly, we can write the magnitude of the isopycnal buoyancy flux ∇b p ·t averaged across surfaces of constant z * as (2.25) The isopycnal buoyancy flux does not appear in the volume-integrated energy budget (2.17)-(2.20) and clearly vanishes for κ S = κ T or when ∇T and ∇S point in the same direction.
In a double-diffusive fluid, an up-gradient buoyancy flux can lead to φ d z * < 0. In practice, it is easier to accurately measure gradients in temperature than salinity in the ocean (Klymak & Nash 2009). It is therefore convenient to recast this condition in terms of the temperature gradient, (2.26) by noting that f (G ρ , θ, 1) = |∇b| 2 /(g 2 β 2 |∇S| 2 ) and applying (2.24). Although G ρ and θ depend on ∇S, if these variables were parametrized, then equation (2.26) could be used with measurements of the temperature gradient to infer the sign of the average diapycnal buoyancy flux and provide insight into the exchange between APE and BPE.
The criteria in equations (2.9) and (2.26) can be applied to turbulent flows, including those in a doubly stable background stratification. We might expect vigorous turbulence to highly distort the scalar surfaces and influence the distribution of θ. If more points lie outside of the range of critical angles such that θ < θ c or θ > 2π − θ c , where θ c is defined in (2.11), then we might anticipate that the average diapycnal buoyancy flux will be positive (down-gradient). This hypothesis will be tested below using simulations of salt fingering.

Evolution of sorted buoyancy
Following the derivation in Winters & D'Asaro (1996) and Nakamura (1996) in the case of a single-component fluid, we can derive an equation for the evolution of the sorted buoyancy profile, averaged in z * coordinates for a double-diffusive fluid: As in the single-component case, the divergence of the diapycnal buoyancy flux sets the time rate of change of the sorted buoyancy. A similar equation was derived by Paparella & von Hardenberg (2012) for the evolution of the horizontally averaged buoyancy profile of a double-diffusive fluid, which evolved due to gradients in a combined advective and diffusive buoyancy flux. One advantage of equation (2.27) is that the diapycnal buoyancy flux does not directly involve the velocity field since the flux is purely diffusive. The single-component version of equation (2.27) was used by Salehipour et al. (2016) to propose a parametrization for the mixing efficiency in z * coordinates, and Taylor & Zhou (2017) used this framework to develop a criterion for layer formation within a single-component stratified flow. Similarly, equation (2.27) could provide a pathway to parametrize double-diffusive flows. For example, it might be possible to parametrize φ d z * as a function of θ and G ρ , although this is beyond the scope of this paper.
As discussed in Nakamura (1996) and Winters & D'Asaro (1996), equation (2.27) can be written in the form of a diffusion equation by defining an effective diffusivity for the sorted buoyancy, κ b (x, t) ≡ φ d z * /(db/dz * ). Therefore, a negative mean diapycnal buoyancy flux implies a negative effective diffusivity, since buoyancy increases monotonically with height in sorted coordinates and hence db/dz * 0. Previous work has shown that double-diffusive flows often develop a negative vertical turbulent diffusivity (e.g. Veronis 1965;St. Laurent & Schmitt 1999;Ruddick & Kerr 2003) and the expression for κ b along with the definition of φ b in equation (2.24) can be viewed as a generalization of this result when the averaging is applied in sorted buoyancy coordinates.

Application: Simulations of salt fingering
To demonstrate an application of the theory described above, we will consider a numerical simulation of a canonical double-diffusive flow: salt fingering. We time-stepped equations (2.1)-(2.5) in a two-dimensional (x-z) domain using a third-order Runge-Kutta scheme and second-order finite differences with an implicit Crank-Nicolson method for the viscous and diffusive terms. Details of the numerical method can be found in Bewley (2012). The initial condition consists of a horizontally uniform temperature and salinity field with linear dependence on height, z. Here the fields are non-dimensionalized such that, for z ∈ (0, 1), the initial temperature field T 0 ∈ (0, 1). Periodic boundary conditions are applied in x with a non-dimensional domain width of 6 to allow multiple salt fingers to form and interact. At the top and bottom of the domain u = 0, and temperature and salinity are set to their initial values. The profiles were set so that warm salty water overlies cold fresh water, a condition favourable for the formation of salt fingers. We used a density ratio of R ρ = 1.01, close in value to 1 to enable the most unstable mode to develop quickly (Kunze 2003). We used a diffusivity ratio of κ T /κ S = 20, which implies θ c ∼ 115 • , and a Prandtl number Pr = ν/κ T = 1 to give a similar critical angle to the oceanic case (θ c ∼ 100) without requiring a high resolution to resolve the Batchelor scales for the scalars. We used typical oceanic values in the linear equation of state (α, β) = (3.9 × 10 −5 • C −1 , 7.9 × 10 −4 ppt −1 ). The simulation was initiated with a small-amplitude random perturbation to the velocity. The results are qualitatively similar for a variety of parameter choices. Figure 2(a) shows the temperature, salinity and buoyancy fields at two times. We only show half of the domain in x for each snapshot; however, the patterns are qualitatively similar across the domain. At t = 300 (left panels) mature salt fingers are visible, and these fingers have broken down into a nonlinear chaotic flow at t = 350 (right panels). The panel in the upper right shows the sorted buoyancy profile, which exceeds the bounds of the initial profile -a situation possible in a double-diffusive fluid. Figure 2(b) shows ∇b p ·n at both times as well as the z * averages. Between the two times shown, the z * average of the diapycnal buoyancy flux changes sign. This appears to be due to the breakdown in the salt fingers, which otherwise preferentially increases salinity gradients.
A random subset of points are superimposed as a scatterplot of ∇b p ·n(G ρ , θ ) on top of f contours in figure 2(c). The unperturbed initial conditions are denoted by a cross. As the salt fingers develop, local gradients in T and S increase. However, |∇S| increases faster than |∇T| and, as a result, G ρ decreases. At t = 300, when salt fingers are present, the scalar angle is primarily distributed around θ = π, but enough points have G ρ < 1 such that the diapycnal buoyancy flux is negative. At t = 350, after the salt fingers have broken down, the distribution of angles becomes bimodal, with clustering of points near θ = π and θ = 0, 2π. The magnitude of the diapycnal buoyancy flux is also maximum near these angles. At this time, the points with ∇b p · n > 0 near θ = 0, 2π more than compensate for the points near θ = π, resulting in a positive mean diapycnal buoyancy flux.
This simulation illustrates how salt fingers distort the temperature and salinity contours such that |∇S| increases faster than |∇T|, leading to G ρ < 1, and BPE conversion into APE on average. Once the flow becomes fully developed, the T and S gradients align such that θ is close to 0, 2π, and as a result the diapycnal buoyancy flux reverses sign and APE is converted into BPE on average. Although there is much left to explore, this example illustrates how the new criteria involving G ρ and θ can be used to analyse the energetics of a double-diffusive flow.

Discussion and conclusion
Here, we extended the energetic framework of Winters et al. (1995), Winters & D'Asaro (1996) and Nakamura (1996) to a double-diffusive fluid. In this framework, an up-gradient diapycnal buoyancy flux averaged in sorted height coordinates is equivalent to a negative buoyancy diffusivity for the sorted buoyancy profile and a transfer of energy from the background potential energy (BPE) to the available potential energy (APE).
We derived criteria for a transfer from BPE to APE in terms of the gradient ratio G ρ = α|∇T|/(β|∇S|), the angle θ between ∇T and ∇S, and the ratio of molecular diffusivities, κ T /κ S (see (2.9) and (2.26)). We applied these critera to salt fingering, an important mixing mechanism within the ocean. The criteria could be applied to ocean observations or direct numerical simulations in other contexts. Finally, we derived an evolution equation for the sorted buoyancy profile.
A transfer of energy from BPE to APE is not possible in a single-component fluid with a linear equation of state, where mixing is an 'irreversible' process. The finding that BPE can be converted into APE within double diffusion implies that some of the BPE can be made 'available' to drive fluid motion, and hence the interpretation of the APE from Winters et al. (1995) does not hold for a double-diffusive fluid. In particular this complicates the definition of 'mixing efficiency' as it is commonly calculated (Gregg et al. 2018). Nevertheless, the criteria in equation (2.26) provides a useful way to quantify the diffusive conversion of potential energy into a form that can be used to drive fluid motion in a double-diffusive fluid.
We formulated the results in this paper assuming molecular diffusion for temperature and salinity. However, a similar analysis could be applied using turbulent diffusivities, i.e. using the equations DT Dt = ∇ · (K T ∇T), DS Dt = ∇ · (K S ∇S), (3.1a,b) where K T (x, t) and K S (x, t) are parametrized turbulent diffusivities. The condition (2.9) for ∇b p ·n < 0 remains the same, with a replacement of molecular diffusivities (κ T , κ S ) with turbulent diffusivities (K T , K S ). Equations (2.17)-(2.20) are unchanged with ∇b p = gαK T ∇T − gβK S ∇S, except that Φ i can no longer be written as a surface term. 893 R3-10 This allows our criterion in equation (2.9) to be applied to ocean models with parametrizations for K T and K S , provided that a linear equation of state remains a good approximation. Although we have focused on fluids where the density is a function of two components, the results hold for an arbitrary buoyancy flux potential, b p . In particular, this implies that our analysis could be extended to include more than two components.