Diapycnal Velocity in the Double-Diffusive Thermocline

A series of large-scale numerical simulations is presented, which incorporate parameterizations of vertical mixing of temperature and salinity by double-diffusion and by small-scale turbulence. These simulations reveal the tendency of double-diffusion to constrain diapycnal volume transport, both upward and downward. For comparable values of mixing coefficients, the average diapycnal velocity in the double-diffusive thermocline is much less than in the corresponding turbulent regime. The insulating effect of double-diffusion is rationalized using two theoretical models. The first argument is based on the assumed vertical advective-diffusive balance. The second theory uses the Rhines and Young technique to evaluate the net diapycnal transport across regions bounded by closed streamlines at a given density surface. The numerical simulations and associated analytical arguments in this study underscore fundamental differences between double-diffusive mixing and mechanically generated small-scale turbulence. When both double-diffusion and turbulence are taken into account, we find that the constraints on diapycnal velocity loosen (tighten) with the increase (decrease) of the fraction of the overall mixing attributed to turbulence. The range of diapycnal velocities that could be realized in doubly-diffusive fluids is determined by the variation in the heat/salt flux ratio. We hypothesize that the unique ability of double-diffusive mixing to actively control diapycnal volume transport may have significant ramifications for the structure and dynamics of thermohaline circulation in the ocean.


Introduction
Theory of heat and mass transfer by small-scale transient motions constitutes a large and important component of numerous physical sciences, including fields as diverse and seemingly unrelated as oceanography, astrophysics, geology and materials science.The key challenge in small-scale mixing theory, in all its forms and applications, is the analysis and prediction of mean turbulent transport for a given large-scale distribution of properties.Most small-scale mixing processes in fluids that are nominally stably stratified-with density decreasing upwards-can be broadly classified into mechanically-driven and double-diffusive categories.An example of a common mechanically-driven phenomenon is the intermittent generation of small-scale turbulence by shear associated with internal gravity waves.Double-diffusive phenomena, on the other hand, are caused by unequal molecular diffusion rates of two (or more) individual density components, such as temperature and salt in the oceanic context.
Recent years have witnessed a substantial progress in our understanding of both double-diffusive and turbulent phenomena.In particular, reliable parameterizations have been developed for double-diffusion in astrophysical [1,2] and oceanographic [3] contexts.These parameterizations express the vertical fluxes (F T , F S ) of diffusing properties (T, S) as a function of the corresponding large-scale gradients using Fick's diffusion model: where the eddy diffusivities (K T , K S ) themselves are not uniform but assumed to depend on the density ratio R ρ = αT z βS z , where (α, β) are the expansion/contraction coefficients.
The purpose of this communication is to discuss one fundamental yet mostly overlooked distinction in the dynamics and consequences of double-diffusive and mechanically-driven mixing.In addition to the vertical transport of diffusing properties (e.g., temperature and salinity) which the flux-gradient laws (Equation ( 1)) attempt to represent, numerous applications require knowledge of the corresponding volume or mass transport across mean density surfaces.Following oceanographic nomenclature, the velocity component normal to the density surfaces (w * ) will be referred to as the diapycnal velocity and the associated volume flux as the diapycnal transport.It is the analysis of diapycnal transport which reveals some rather unique properties of double-diffusion that have no direct counterpart in one-component turbulent flows.In particular, we argue that in double-diffusive flows, diapycnal volume transport is dramatically reduced relative to that in turbulent systems with comparable mixing rates of density components.When both double-diffusion and turbulence are present, diapycnal velocity is restricted to a finite range, which is controlled by the relative contributions of double-diffusion and turbulence to overall mixing.
Much of the motivation for the analysis of diapycnal volume transport [4][5][6][7] stems from the role it plays in the theory of thermohaline circulation in the ocean.Thermohaline circulation-the system of large-scale flows driven by the sea-water density variation-is an essential regulator of the Earth's climate.Nevertheless, several aspects of it have not been fully explained.The classical view [8][9][10][11] emphasizes the role of small-scale mixing and associated water mass transformation in the ocean interior.This idea was challenged by studies advocating an alternative adiabatic view [12][13][14] which invokes isopycnal advection as the key mechanism for maintenance of the Meridional Overturning Circulation (MOC).While the exact contributions of different circulation components remain a source of debate, there is a general consensus that the diabatic mode, driven by diapycnal fluxes, can be substantial.Radko et al. [15] for instance, estimate that the diabatic mode of circulation accounts for approximately one third of the Meridional Overturning Circulation (MOC) in the global ocean.Another uncertain aspect of the theory of thermohaline circulation concerns the relative role of shallow overturning cells in the main thermocline and deep circulation in the abyssal ocean.Earlier efforts [9,16] were mostly focused on the abyssal dynamics.However, an interesting suggestion was made [17] that the meridional transport of heat is controlled by the processes operating in the upper ocean.In an attempt to quantify partitioning of the oceanic heat flux between the deep and shallow overturning cells, these authors introduced the "heatfunction", a quantity which identifies the components of circulation which contribute to the total meridional heat flux.Boccaletti et al. [17] have shown that the heatfunction represents a surface-intensified flow largely limited to the main thermocline.Therefore meridional heat transport is dominated by the contribution from the upper branch of the MOC where double-diffusion is most active.For instance, as much as 95% of the Atlantic thermocline is double-diffusively unstable [18].Thus, the link between double-diffusion and diapycnal volume transport can have potentially critical ramifications for thermohaline circulation and should be explored much further.
In order to rationalize the tendency of double-diffusion to suppress the diapycnal volume transport, we propose two arguments.The first argument applies to the strongly diabatic regime representative of the lower (diffusive) thermocline that is often described by Munk's [9] vertical advection-diffusion balance.We find that steady one-dimensional solutions in double-diffusive systems are possible only for a very narrow range of diapycnal velocities.The second argument is more relevant for the central thermocline, where double-diffusive processes are most common.Dynamics of the main thermocline are fundamentally three-dimensional [19].The overall water-mass distribution is set by large-scale advection whereas diapycnal mixing is relatively weak.Nevertheless, in this regime it is still possible to predict the impact of higher order diabatic processes on cross-isopycnal transfer using a technique originally developed by Rhines and Young [20].The crux of this approach is the integration of governing equations along closed streamlines, which effectively eliminates the zero order advective terms in integral balances.Adapting this procedure to the double-diffusive problem allows us to evaluate the average diapycnal velocity in regions bounded by closed streamlines and, ultimately, to explain the insulation effect from first principles.Thus, for both advectively-dominated and diffusively-dominated regimes, there are theoretical reasons to expect the suppression of diapycnal volume flux by double-diffusion.
The paper is set as follows.Our starting point (Section 2) is an analysis of a one-dimensional vertical advective-diffusive balance of temperature and salinity of Munk's [9] type.This analysis, albeit highly idealized, suggests that fundamental differences may exist in the way turbulent and double-diffusive mixing affect diapycnal volume transport.Motivated by this possibility, we perform a series of large-scale three-dimensional multi-century numerical simulations (Section 3), which also consistently reflect the tendency of double-diffusion to constrain diapycnal volume transport.In Section 4, we develop an analytical model of the double-diffusive insulation in a three-dimensional setting.We summarize and conclude in Section 5.

Preliminary Considerations
Our inquiry into the link between double-diffusive mixing and the associated diapycnal volume transport was originally motivated by the following fairly abstract problem.Consider a one-dimensional model consisting of the steady state temperature and salinity equations: where T(z) and S(z) are the large-scale temperature and salinity of the sea-water; w is the upward vertical velocity.The temperature and salinity fluxes (F T , F S ) are attributed to small-scale mixing processes and related to large-scale gradients through Equation (1).We are interested in solving Equation (2) for given boundary conditions: where z = −H top , (z = −H bot ) represents the top (bottom) of the mixing zone and H is its thickness (see the schematic in Figure 1.In particular, we are concerned by the differences in the solutions of Equations ( 2) and (3) for cases in which vertical mixing is controlled by: mechanical turbulence, double-diffusion, and a combination thereof.).

Turbulent Mixing
To illustrate the basic properties of the 1D system (Equation ( 2)), consider first the simplest model of turbulent mixing characterized by constant and equal eddy diffusivities of heat and salt: In this case, Equations ( 2) and (3) can be solved for any value of diapycnal velocity ( w ): and Thus, the turbulent mixing model (Equation ( 4)) represents an extreme case of a completely passive system which offers no restrictions on diapycnal velocity.For any value of w it is possible to construct a consistent solution of the T-S advection-diffusion equations satisfying the boundary conditions at ).The latter model can be used to represent differential diffusion, the preferential turbulent mixing of temperature relative to salinity under certain conditions [21,22].In order to demonstrate that the passive character of the advective-diffusive balance is not necessarily realized in other types of mixing, we now turn to a very different configuration which demands a unique value of w-the constant flux ratio model of double-diffusion.

Double-Diffusion: The Constant Flux Ratio Model
Double-diffusion comes in two distinct forms, fingering and diffusive convection.The fingering regime is commonly realized in the subtropical thermocline where warm and salty water overlies cold and fresh.Diffusive convection is more prevalent in high-latitude oceans, where temperature

Turbulent Mixing
To illustrate the basic properties of the 1D system (Equation ( 2)), consider first the simplest model of turbulent mixing characterized by constant and equal eddy diffusivities of heat and salt: In this case, Equations ( 2) and (3) can be solved for any value of diapycnal velocity (w): Thus, the turbulent mixing model (Equation ( 4)) represents an extreme case of a completely passive system which offers no restrictions on diapycnal velocity.For any value of w it is possible to construct a consistent solution of the T-S advection-diffusion equations satisfying the boundary conditions at z = z bot and z = z top .It should be noted that the advective-diffusive balance (Equation ( 2)) remains passive even when the diffusivities of heat and salt are spatially uniform but not equal to each other (K T = K S ).The latter model can be used to represent differential diffusion, the preferential turbulent mixing of temperature relative to salinity under certain conditions [21,22].In order to demonstrate that the passive character of the advective-diffusive balance is not necessarily realized in other types of mixing, we now turn to a very different configuration which demands a unique value of w-the constant flux ratio model of double-diffusion.

Double-Diffusion: The Constant Flux Ratio Model
Double-diffusion comes in two distinct forms, fingering and diffusive convection.The fingering regime is commonly realized in the subtropical thermocline where warm and salty water overlies cold and fresh.Diffusive convection is more prevalent in high-latitude oceans, where temperature and salinity frequently increase downward.While we believe that the analysis in this paper is relevant for both forms of double-diffusion, to be specific, the discussion is focused on the salt-finger regime.
The single most significant parameter controlling diapycnal mixing in regions of active double-diffusion is the local density ratio where α and β are the expansion/contraction coefficients in the linear equation of state-see the discussion in Schmitt [18,23] and Radko [24].Double-diffusive parameterizations typically assume that the eddy diffusivities of heat and salt (K T , K S ) are determined by R ρ .It is also common to approximate the ratio of thermal and haline density fluxes, so-called flux ratio (γ), by a constant value: Particular examples of such parameterizations can be found in, for example, Schmitt [25], Zhang et al. [26], and Merryfield et al. [27].
Combining Equations ( 2) and ( 8), we arrive at To maintain fingering convection, the loss of potential energy stored in salt stratification should exceed the energy gain by heat stratification and therefore γ 0 < 1.The density ratio, on the other hand, has to exceed unity in order for density to decrease upward.Hence, R ρ γ 0 > 1 and therefore Equation ( 9) can be satisfied only if w = 0. Thus, in contrast with the turbulent mixing model, the double-diffusive system actively controls diapycnal velocity.Of course, the selection of a unique w in this model could be a consequence of the chosen flux laws; of particular concern is the constant flux ratio approximation (Equation ( 8)).However, even when the assumption of constant flux ratio is relaxed (Appendix) by considering typical finite variations in γ(R ρ ) we still find that the range of permitted diapycnal velocities is rather narrow: max |w| ∼ 5 × 10 −9 m/s.

Combined Effects of Double-Diffusion and Turbulence
Because small-scale mixing in the ocean is controlled by a combination of double-diffusion and mechanically generated turbulence, it is of interest to examine the constraints on diapycnal velocity in a model which includes both mixing processes.For that, we extend the foregoing analyses by considering the following closure: which assumes equal and uniform diffusivities of heat and salt due to turbulence (K turb ).Note that this simplified conceptual model does not take into account the possibility of direct influence of turbulence on double-diffusion [28,29].
There have been at least two attempts to evaluate the flux ratio based on oceanographic field measurements [30,31].Both estimates were mutually consistent, suggesting the flux ratio of γ ≈ 0.85 (11) The individual values of salt finger fluxes, and particularly the pattern of their variation with R ρ , are more difficult to infer from observations.Therefore, we adopt the parameterization proposed by Radko and Smith [3] on the basis of high-resolution DNS: where (a S , b S , R cuto f f ) = (135.7,−62.75, 5.67) and We assume parameter values representative of mid-latitude thermocline: where R 0 = α∆T β∆S .Next, we systematically vary K turb in an attempt to determine how it affects the range of diapycnal velocities.
The system consisting of Equation ( 2) and Equations ( 11)-( 14) was coded in Maple software package and solutions were sought for a range of w, K turb .Note that in the context of the one-dimensional model, in which all field variables are assumed to vary in z only, the continuity equation ∇ • → v = 0 reduces to the requirement that w = const.This feature simplifies the exploration of the parameter space and the results are shown in Figure 2. Calculations resulting in regular solutions are indicated by the heavy dots and the light dots represent conditions under which no solutions were found.Figure 2 indicates that the constraints on diapycnal velocity loosen (tighten) with the increase (decrease) in K turb /K dd T (R 0 )-a quantity measuring the relative contributions from turbulence and double-diffusion to the net mixing.The individual values of salt finger fluxes, and particularly the pattern of their variation with R ρ , are more difficult to infer from observations.Therefore, we adopt the parameterization proposed by Radko and Smith [3] on the basis of high-resolution DNS: where ( , , ) (135.7, 62.75, 5.67) We assume parameter values representative of mid-latitude thermocline: where . Next, we systematically vary turb K in an attempt to determine how it affects the range of diapycnal velocities.
The system consisting of Equation ( 2) and Equations ( 11)-( 14) was coded in Maple software package and solutions were sought for a range of ( ) , turb w K . Note that in the context of the one-dimensional model, in which all field variables are assumed to vary in z only, the continuity equation v 0 ∇⋅ =  reduces to the requirement that w c o n s t = .This feature simplifies the exploration of the parameter space and the results are shown in Figure 2. Calculations resulting in regular solutions are indicated by the heavy dots and the light dots represent conditions under which no solutions were found.Figure 2 indicates that the constraints on diapycnal velocity loosen (tighten) with the increase (decrease) in    Despite considerable uncertainties in the estimated magnitudes of double-diffusive and turbulent mixing, it is plausible that in the central thermocline of the Atlantic Ocean their diffusivities are comparable [32].For K turb /K dd T (R 0 ) ∼ 1 we obtain the following constraint on the (dimensional) diapycnal velocity: max |w| ∼ 10 −7 m/s ∼ 3 m/year (15) which is similar to the data-based estimates of the diapycnal velocity in regions susceptible to both double-diffusive and turbulent mixing [32].In summary, the foregoing calculations indicate that (i) in the turbulent model, the assumed vertical T-S balances do not impose any internal constraints on diapycnal transport; (ii) in the double-diffusive model the diapycnal transport is zero; and (iii) the hybrid model, which takes into account both double-diffusion and turbulence, allows only a finite range of diapycnal velocities.These findings, of course, should be interpreted with great caution-the ability of any one-dimensional model to reflect the inherently three-dimensional dynamics of ocean circulation is suspect.Perhaps it is most profitable to interpret these one-dimensional solutions as the result of the averaging of the advection-diffusion equations along isopycnal surfaces, although this conceptualization raises the question whether such averaging is physically justified.Nevertheless, the profound differences in the way various mixing models affect diapycnal volume transport provide a compelling reason to launch a more comprehensive investigation of this phenomenon using more realistic-three-dimensional and time-dependent-models.We now proceed with the analysis of double-diffusive insulation based on a series of large-scale numerical simulations.

Formulation
The MIT General Circulation Model (MITgcm), described in [33,34], was used to simulate an ocean of equatorial to sub-arctic meridional extent with a zonal width comparable to that of the North Atlantic.
where τ 0 = 0.05 N • m −2 .The variation in the Coriolis parameter (f ) is represented by the beta-plane approximation: where The model was initialized from rest with a T-S distribution that is zonally uniform but varies linearly in y and z.The target temperatures and salinities (Figure 3) served as initial surface values as well and initial bottom values were ⋅ .The nonlocal K-Profile Parameterization (KPP) scheme of Large et al. [35] is used to model convective and near-surface mixing processes.Geostrophic eddy parameterization follows Gent and McWilliams [36].To reduce spurious interactions with KPP and to isolate the dynamics of interest, the isopycnal diffusivity was set to the minimal value ( Vertical mixing of temperature and salinity was represented by a combination of double-diffusion and turbulence.Double-diffusion was parameterized using Equations ( 12) and ( 13).The interior turbulent diffusivity ( turb K ) associated with overturning gravity waves, is assumed to be spatially uniform and equal for temperature and salinity.In a series of simulations ⋅ .This restriction was implemented to prevent excessive convection in the regions where the western boundary current, transporting relatively warm and light water northward, enters locations with much higher target surface density.Each simulation was integrated forward for at least 200 years using a time step of 480 s, which was sufficient to produce quasi-steady circulation patterns.
Figures 4 and 5 present a typical final state realized in our numerical experiments.This simulation has been performed assuming purely double-diffusive interior mixing, with Figure 4a shows the horizontal temperature distribution and the velocity pattern at 500 m z = − .As expected, the flow field is dominated by clockwise circulation in the subtropical region ( ) where wind forcing is anticyclonic.A vertical zonal section of temperature at 0.5 y y L = is shown in Figure 4b.Clearly visible is a well-defined thermocline with relatively warm water extending several hundred meters downward from the surface.and max ρ ( min ρ ) represents the maximal (minimal) density in the computational domain.The color coding represents the values of the density ratio on this isopycnal.The observed range of density ratios < indicates that fingering is active on this isopycnal and that its intensity is relatively high.[36].To reduce spurious interactions with KPP and to isolate the dynamics of interest, the isopycnal diffusivity was set to the minimal value (κ ρ = 10 m 2 • s −1 ) required to maintain the numerical stability of the model.Vertical mixing of temperature and salinity was represented by a combination of double-diffusion and turbulence.Double-diffusion was parameterized using Equations ( 12) and (13).The interior turbulent diffusivity (K turb ) associated with overturning gravity waves, is assumed to be spatially uniform and equal for temperature and salinity.In a series of simulations K turb was systematically varied from 0 to 2 × 10 −4 m 2 • s −1 .In the region of western intensification (x < L int ), however, K turb was maintained at the minimum value of K turb min WBC = 2 × 10 −4 m 2 • s −1 .This restriction was implemented to prevent excessive convection in the regions where the western boundary current, transporting relatively warm and light water northward, enters locations with much higher target surface density.Each simulation was integrated forward for at least 200 years using a time step of 480 s, which was sufficient to produce quasi-steady circulation patterns.
Figures 4 and 5 present a typical final state realized in our numerical experiments.This simulation has been performed assuming purely double-diffusive interior mixing, with K turb = 0 for x > L int = 1000 km. Figure 4a shows the horizontal temperature distribution and the velocity pattern at z = −500 m.As expected, the flow field is dominated by clockwise circulation in the subtropical region ( 1 / 4 L y < y < 3 / 4 L y ) where wind forcing is anticyclonic.A vertical zonal section of temperature at y = 0.5L y is shown in Figure 4b.Clearly visible is a well-defined thermocline with relatively warm water extending several hundred meters downward from the surface.Figure 5 presents a three-dimensional view of the isopycnal ρ = ρ av , where and ρ max (ρ min ) represents the maximal (minimal) density in the computational domain.The color coding represents the values of the density ratio on this isopycnal.The observed range of density ratios 1 < R ρ < 3 indicates that fingering is active on this isopycnal and that its intensity is relatively high.Color coding represents the density ratio distribution on this isopycnal.The observed range of density ratios 1 3 R ρ < < indicates that fingering is active and that its intensity is relatively high.

Diapycnal Velocity
The analysis in Section 2 pertained to steady one-dimensional systems in which diapycnal volume transport is associated with vertical advection.To quantify diapycnal fluxes in a more realistic three-dimensional setting, the diagnostics of large-scale simulations will be focused on diapycnal velocity * w [37,38] which is generally dissimilar to the vertical velocity ( w).Diapycnal velocity can be defined by rewriting governing equations in terms of density (rather than z) as a vertical coordinate: where In Equation ( 21),     Color coding represents the density ratio distribution on this isopycnal.The observed range of density ratios 1 3 R ρ < < indicates that fingering is active and that its intensity is relatively high.

Diapycnal Velocity
The analysis in Section 2 pertained to steady one-dimensional systems in which diapycnal volume transport is associated with vertical advection.To quantify diapycnal fluxes in a more realistic three-dimensional setting, the diagnostics of large-scale simulations will be focused on diapycnal velocity * w [37,38] which is generally dissimilar to the vertical velocity ( w).Diapycnal velocity can be defined by rewriting governing equations in terms of density (rather than z) as a vertical coordinate: where In Equation ( 21), ( )  Color coding represents the density ratio distribution on this isopycnal.The observed range of density ratios 1 < R ρ < 3 indicates that fingering is active and that its intensity is relatively high.

Diapycnal Velocity
The analysis in Section 2 pertained to steady one-dimensional systems in which diapycnal volume transport is associated with vertical advection.To quantify diapycnal fluxes in a more realistic three-dimensional setting, the diagnostics of large-scale simulations will be focused on diapycnal velocity w * [37,38] which is generally dissimilar to the vertical velocity (w).Diapycnal velocity can be defined by rewriting governing equations in terms of density (rather than z) as a vertical coordinate: where In Equation (21) Figure 6 illustrates the following physical interpretation of diapycnal velocity.Consider a Lagrangian particle initially located at the motionless isopycnal surface ρ (point A).In the presence of a finite cross-isopycnal flow, this particle will be advected by velocity → v = (u, v, w) to a new point B, located off the isopycnal.Let C be the vertical projection of the point B onto the isopycnal surface.In this configuration, the diapycnal velocity in Equation (22) represents the rate of increase of the vertical separation of our Lagrangian point from the density surface: Note that the shadow point C, which remains on the isopycnal, can also alter its z-level and the corresponding vertical velocity can be expressed as Thus, the isopycnal ascent rate w ρ and diapycnal velocity w * represent two distinct components of vertical velocity: Figure 6 illustrates the following physical interpretation of diapycnal velocity.Consider a Lagrangian particle initially located at the motionless isopycnal surface ρ (point A).In the presence of a finite cross-isopycnal flow, this particle will be advected by velocity v ( ,v, ) u w =  to a new point B, located off the isopycnal.Let C be the vertical projection of the point B onto the isopycnal surface.In this configuration, the diapycnal velocity in Equation (22) represents the rate of increase of the vertical separation of our Lagrangian point from the density surface: Equation (22) indicates that in the absence of any diabatic water-mass transformation, and therefore diapycnal velocity is zero.Curiously, the reverse statement is not necessarily correct.The double-diffusive example in Section 2.2 indicates that diapycnal velocity can still be zero even in the presence of substantial diabatic mixing of properties.An alternative explanation of low diapycnal velocity can be based on the dianeutral velocity equation [4] which, for the conditions assumed in our study, reduces (in present notation) to: In the case of a purely double-diffusive system, the diapycnal transport is represented by the third term on the right-hand side, which is proportional to ( ) Since the flux ratio in the ocean (Equation ( 11)) is rather close to unity, relatively low values of diapycnal velocity are expected even when the stratification cannot be accurately described by the predominantly one-dimensional Munk-type model.Equation (22) indicates that in the absence of any diabatic water-mass transformation, dρ dt = → v • ∇ρ = 0 and therefore diapycnal velocity is zero.Curiously, the reverse statement is not necessarily correct.The double-diffusive example in Section 2.2 indicates that diapycnal velocity can still be zero even in the presence of substantial diabatic mixing of properties.
An alternative explanation of low diapycnal velocity can be based on the dianeutral velocity equation [4] which, for the conditions assumed in our study, reduces (in present notation) to: In the case of a purely double-diffusive system, the diapycnal transport is represented by the third term on the right-hand side, which is proportional to (1 − γ).Since the flux ratio in the ocean (Equation ( 11)) is rather close to unity, relatively low values of diapycnal velocity are expected even when the stratification cannot be accurately described by the predominantly one-dimensional Munk-type model.
The following analysis will therefore attempt to determine whether, and to what extent, w * is suppressed by double-diffusion in the more realistic three-dimensional systems.

Diapycnal Transport in the Double-Diffusive and in the Turbulent Ocean
First, the proposed diagnostics have been applied to the double-diffusive experiment shown in Figures 4 and 5. Diapycnal velocity was computed using Equation ( 22) and evaluated at the isopycnal surface ρ = ρ av .The results are plotted (Figure 7a) for the interior region (L int < x < L x ) where mixing is exclusively double-diffusive.The same procedure was then applied to the turbulent experiment (Figure 7b).The latter was designed in the same way as the baseline double-diffusive run, except that the vertical mixing was assumed to be exclusively turbulent with The comparison of diapycnal velocities in these experiments reflects the dramatic difference in the way double-diffusion and turbulence affect diapycnal transport.Over most of the isopycnal surface, the double-diffusive w * are on the order of ∼ 5 × 10 −9 m/s or less (both positive and negative values are observed).On the other hand, the turbulent w * is mostly positive and is larger by at least an order of magnitude (∼ 5 × 10 −8 m/s).The difference between Figure 7a  The following analysis will therefore attempt to determine whether, and to what extent, * w is suppressed by double-diffusion in the more realistic three-dimensional systems.

Diapycnal Transport in the Double-Diffusive and in the Turbulent Ocean
First, the proposed diagnostics have been applied to the double-diffusive experiment shown in Figures 4 and 5. Diapycnal velocity was computed using Equation ( 22) and evaluated at the isopycnal surface av ρ ρ = .The results are plotted (Figure 7a) for the interior region ( int x L x L < < ) where mixing is exclusively double-diffusive.The same procedure was then applied to the turbulent experiment (Figure 7b).The latter was designed in the same way as the baseline double-diffusive run, except that the vertical mixing was assumed to be exclusively turbulent with  Figure 8a presents the mean diapycnal velocity evaluated at various isopycnal surfaces w * iso (ρ).In computing w * iso , we excluded the contribution from two regions with elevated turbulent diffusivity: (i) the western intensification region (x < L int ) and (ii) the near surface region (z > −250 m) where some isopycnal surfaces enter the mixed layer, particularly in the cold Northern parts of the domain.Diapycnal velocity w * iso was evaluated for both turbulent and double-diffusive experiments and referenced (Figure 8a) to the average depth of isopycnal surfaces h iso (ρ). Figure 8a shows that turbulent diapycnal velocity (blue curve) substantially exceeds the double-diffusive w * iso (green curve) for most isopycnals.In both cases, relatively large values of w * iso were found at the base of the main thermocline (h av ∼ 1000 m).Diapycnal velocity in the turbulent experiment was also greatly amplified in the abyssal region ( h → 3000 m ) whereas no such amplification occurred in the double-diffusive simulation.It is interesting to note that the patterns of isopycnal-averaged diapycnal velocity are qualitatively similar to the corresponding local vertical profiles.For instance, in Figure 8b we plot the diapycnal velocity profile at the center of the model domain: w * loc = w * (0.5L x , 0.5L y , z).Local diagnostics also indicate that the diapycnal velocity in the turbulent ocean greatly exceeds the double-diffusive w * .Thus, the double-diffusive insulation effect considered in this study represents both the integral and the local property of large-scale flows.
double-diffusive simulation.It is interesting to note that the patterns of isopycnal-averaged diapycnal velocity are qualitatively similar to the corresponding local vertical profiles.For instance, in Figure 8b we plot the diapycnal velocity profile at the center of the model domain: An attempt has also been made to determine the constraints on diapycnal velocity in an ocean experiencing vertical mixing due to a combination of double-diffusion and mechanically generated turbulence.For that, we performed a series of simulations in which turb K was systematically varied, while retaining the parameterization Equations ( 11)-( 13) for double-diffusion.Diapycnal velocity was computed and evaluated as previously (Figures 7 and 8).As expected, we find that diapycnal transport intensifies with the increase of the fraction of the overall mixing attributed to turbulence.An attempt has also been made to determine the constraints on diapycnal velocity in an ocean experiencing vertical mixing due to a combination of double-diffusion and mechanically generated turbulence.For that, we performed a series of simulations in which K turb was systematically varied, while retaining the parameterization (Equations ( 11)-( 13)) for double-diffusion.Diapycnal velocity was computed and evaluated as previously (Figures 7 and 8).As expected, we find that diapycnal transport intensifies with the increase of the fraction of the overall mixing attributed to turbulence. Figure 9 presents the variation of the isopycnal-averaged diapycnal velocity w * av = w * iso (ρ av ) at the isopycnal surfaces ρ = ρ av in Equation ( 19) as a function of K turb .These diagnostics also support our interpretation of double-diffusive mixing as an insulating mechanism.When K turb is less than typical values for K dd S and K dd T , the overall system behaves like the fully double-diffusive case (e.g., Figure 7a).Values of mean diapycnal velocity are on the order of ∼ 5 × 10 −9 m/s or less for experiments with K turb ≤ 2 × 10 −6 m 2 /s.Conversely, simulations with K turb ≥ K dd S , K dd T ∼ 2 × 10 −5 m 2 /s are characterized by w * on the order of ∼ 5 × 10 −8 m/s or greater.

The Role of the Flux Ratio
The foregoing experiments have illustrated the dramatic differences in diapycnal volume flux in double-diffusive and turbulent systems.The question that naturally arises at this point is what aspect of double-diffusion makes it so effective in blocking diapycnal volume transport?Theoretical

The Role of the Flux Ratio
The foregoing experiments have illustrated the dramatic differences in diapycnal volume flux in double-diffusive and turbulent systems.The question that naturally arises at this point is what aspect of double-diffusion makes it so effective in blocking diapycnal volume transport?Theoretical arguments presented in the Appendix suggest that, at least in one-dimensional models, the range of diapycnal velocities is controlled by the effective variation in the flux ratio (γ).Double-diffusive systems are characterized by a relatively uniform γ.For instance, Radko and Smith [3] find that as the density ratio increases from R ρ = 1.1 to R ρ = 3, the flux ratio decreases by ∆γ ∼ 0.15.The flux ratio realized in turbulent fluids γ = R ρ varies over the same interval much more (∆γ = 1.9), which suggests significantly larger values of diapycnal velocity.In particular, the one-dimensional model (Appendix) predicts that w * ∝ ∆γ.To determine whether these conclusions remain relevant for three-dimensional systems, we have performed an additional series of four double-diffusive simulations with variable flux ratio patterns.
The new simulations were identical to the baseline double-diffusive experiment (Figures 3 and 4) in all respects except for the chosen flux ratio model, which assumed that γ varies linearly with the density ratio: γ The flux ratio models used in these simulations are shown in Figure 10a, along with the baseline (δ = 0) experiment, and the corresponding parameter values are δ = 0, 0.25, 0.5, 0.75, 1 and γ 0 = 0.85, 0.8875, 0.925, 0.9625, 1 respectively.The negative values of δ have not been considered since the decreasing γ R ρ relation results in the formation of thermohaline staircases [39][40][41].The presence of staircases, characterized by nearly homogeneous interior convecting layers, would greatly complicate the diagnostics of diapycnal velocity (Equation ( 22)).Thus, the flux ratio pattern in Figure 10a gradually changes from the zero-slope baseline experiment (γ = 0.85) to the one realized in the turbulent model γ = R ρ .In all experiments, however, we retain the double-diffusive mixing parameterization (Equations ( 12) and ( 13)); the corresponding K S (R ρ ) pattern is shown in Figure 10b.As previously (e.g., Figure 9), our analysis is focused on the properties of * broadly consistent with the one-dimensional theoretical model developed in the Appendix (see Figure A2).As previously (e.g., Figure 9), our analysis is focused on the properties of w * av -the isopycnal-averaged diapycnal velocity calculated for each simulation along the density surface ρ = ρ av -and these results are summarized in Figure 11.The diapycnal velocity for the δ = 0 case is very weak (w * ∼ 10 −9 m/s) and negative.However, as δ increases and the flux ratio becomes less and less uniform, w * av monotonically increases.Finally, for δ = 1, which corresponds to the turbulent flux ratio model γ = R ρ , the average diapycnal velocity increases to w * av ≈ 5 × 10 −8 m/s.
Diapycnal velocities in this case become comparable to those realized in the fully turbulent model (Figure 7b).The pattern of w * av (δ) in Figure 11 can be reasonably well described as a linear relation and is broadly consistent with the one-dimensional theoretical model developed in the Appendix (see Figure A2).
As previously (e.g., Figure 9), our analysis is focused on the properties of * .Diapycnal velocities in this case become comparable to those realized in the fully turbulent model (Figure 7b).The pattern of * ( δ ) a v w in Figure 11 can be reasonably well described as a linear relation and is broadly consistent with the one-dimensional theoretical model developed in the Appendix (see Figure A2).

Forced Simulations
The previous experiments addressed the evolution of diapycnal transport driven by surface wind and thermohaline forcing through relaxation of temperature and salinity at the surface.To explore the effectiveness of double-diffusion in constraining diapycnal velocity under more controlled conditions, another series of "forced" experiments was performed in which the vertical velocity was set at the target value above the sea-floor and below the ocean surface.The interior Figure 11.A series of large-scale simulations in which the variation in the flux ratio (as measured by the parameter δ) is systematically increased.For each experiment, the mean diapycnal velocity (w * av ) at the average isopycnal surface ρ av is plotted as a function of δ.Note the monotonic-nearly linear-increase in w * av with δ.

Forced Simulations
The previous experiments addressed the evolution of diapycnal transport driven by surface wind and thermohaline forcing through relaxation of temperature and salinity at the surface.To explore the effectiveness of double-diffusion in constraining diapycnal velocity under more controlled conditions, another series of "forced" experiments was performed in which the vertical velocity was set at the target value above the sea-floor and below the ocean surface.The interior upwelling was implemented in MITgcm by relaxing the horizontal velocities in the surface and bottom boundary layers to the prescribed convergent (at the bottom) and divergent (below the surface) patterns.The interior upwelling was balanced by downwelling in the vicinity of the Northern boundary.Wind forcing was not applied.One of the goals of the forced experiments was the determination of system behavior when the target vertical velocity significantly exceeds the anticipated constraint on diapycnal velocity.
The first experiment was performed using the fully turbulent mixing model with K T = K S = 10 −4 m 2 /s-the nominal Munk [9] diffusivity-and the target value for vertical velocity was 10 −7 m/s.The vertical T-S profiles at the center of the model domain 0.5L x , 0.5L y after 200 years of integration are shown in Figure 12a.As expected for the turbulent model, which offers no constraints on the magnitude of diapycnal velocity, stratification takes a regular exponential form.Figure 12b presents the corresponding double-diffusive simulation.In this experiment, we used the functional form (Equation ( 12)) for the parameterization of double-diffusion, but the amplitude was adjusted upward to match the diffusivity assumed for the turbulent experiment (Figure 12a).The outcomes of these two simulations are drastically different-the stratification in the double-diffusive simulation becomes highly irregular, with fingering regions separated by deep well-mixed convecting layers.
It is also interesting to note that values of diapycnal velocity evaluated in the double-diffusive regions remain small-at least an order of magnitude less than in the corresponding turbulent experiment (Figure 12a) and, somewhat surprisingly, significantly less than the target vertical velocity.These results underscore the robust and generic tendency of double-diffusion to constrain diapycnal velocity in various configurations.
Figure 12b presents the corresponding double-diffusive simulation.In this experiment, we used the functional form (Equation ( 12)) for the parameterization of double-diffusion, but the amplitude was adjusted upward to match the diffusivity assumed for the turbulent experiment (Figure 12a).The outcomes of these two simulations are drastically different-the stratification in the double-diffusive simulation becomes highly irregular, with fingering regions separated by deep well-mixed convecting layers.It is also interesting to note that values of diapycnal velocity evaluated in the double-diffusive regions remain small-at least an order of magnitude less than in the corresponding turbulent experiment (Figure 12a) and, somewhat surprisingly, significantly less than the target vertical velocity.These results underscore the robust and generic tendency of double-diffusion to constrain diapycnal velocity in various configurations.

Theoretical Model of Double-Diffusive Insulation
Numerical simulations (Section 3) indicate that double-diffusive insulation occurs over most of the thermocline and it is not limited to Munk's regions characterized by vertical advective-diffusive balance.Therefore our next step represents an attempt to rationalize this phenomenon without relying on the highly restrictive one-dimensional assumption used in Section 2.Over much of the main thermocline, small-scale mixing constitutes a numerically small component of the full advective-diffusive balance of temperature and salinity.Yet, this component is responsible for the diabatic water-mass transformation and, ultimately, for the selection of diapycnal velocity.Diabatic

Theoretical Model of Double-Diffusive Insulation
Numerical simulations (Section 3) indicate that double-diffusive insulation occurs over most of the thermocline and it is not limited to Munk's regions characterized by vertical advective-diffusive balance.Therefore our next step represents an attempt to rationalize this phenomenon without relying on the highly restrictive one-dimensional assumption used in Section 2.Over much of the main thermocline, small-scale mixing constitutes a numerically small component of the full advective-diffusive balance of temperature and salinity.Yet, this component is responsible for the diabatic water-mass transformation and, ultimately, for the selection of diapycnal velocity.Diabatic dynamics can be brought to the fore and analyzed theoretically using a procedure analogous to the one used by Rhines and Young [20], albeit in a very different context.
Consider a zero-order steady state of the ocean in the absence of mixing.This state is characterized by exact conservation of temperature, salinity, density and potential vorticity.This ideal basic state is slightly perturbed by including weak diapycnal fluxes of heat and salt in the advection-diffusion equations of motion.It is assumed that the perturbation results in only slight (first order) modification of the corresponding ideal zero order solution.The perturbed state is governed by where and w ρ is given in Equation (23). Figure 13 represents the typical circulation pattern in the subtropical ocean on an arbitrarily chosen isopycnal surface ρ = ρ 0 .On this surface, we select a closed temperature contour (T = T 0 ).Since temperature, salinity and density are assumed to be uniquely related through the linear equation of state, it follows that salinity is also constant along the chosen contours (S = S 0 ).These contours closely follow the corresponding streamlines of the ideal basic state on which (T 0 , S 0 , ρ 0 ) are conserved exactly.The following analysis is pivoted about the time integrals of Equation ( 28) along such closed contours: The major simplification here is achieved by elimination of the zero order advective terms.Combining Equation (29) with the double-diffusive parameterization (Equation ( 8)), we arrive at The following analysis is pivoted about the time integrals of Equation ( 28) along such closed contours: The major simplification here is achieved by elimination of the zero order advective terms.Combining Equation (29) with the double-diffusive parameterization (Equation ( 8)), we arrive at The next step is the change of the integration variable from time (t) to the length coordinate measured along the contours (l), which reduces Equation (30) to where → v h is the absolute value of horizontal velocity.The mid-latitude circulation patterns (e.g., Figure 12) consist of two distinct regions: relatively slow geostrophic interior (L int < x < L x ) and swift western boundary currents (0 < x < L int ).The inspection of Equation (31) indicates that even if a streamline passes through both regions (Figure 13), the dominant contribution to Equation (31) comes from the interior geostrophic region, where → v h is much less (by a factor L x /L WBC >> 1) than the velocity in the boundary current.Neglecting the small contribution from the region of western intensification allows us to focus exclusively on geostrophic dynamics and express the horizontal velocity components as follows: where M is the Montgomery potential.This system (Equation ( 32)) indicates that the contours of integration in Equation ( 31) at the leading order coincide with the isopleths of M. Another simplification brought by geostrophic approximation is that the Ertel potential vorticity in this case takes a relatively simple form where C is any conservative tracer.Recall that the perturbed solution is assumed to be only slightly different from its adiabatic counterpart, which conserves temperature, salinity, and potential vorticity in the Lagrangian sense.Thus, we conclude that at the leading order (i) potential vorticity q is conserved along the contours of integration and (ii) the tracer C can be represented by any combination of temperature and salinity.Given the structure of the line integral (Equation ( 31)), we expect significant simplifications to occur for Since potential vorticity is approximately uniform along each streamline, we divide Equation (31) by q, entering it directly into the integrand, which reduces the line integral to Here, the notation "int" is used to emphasize that the integration is carried out only over the geostrophic interior.By virtue of Equation ( 32), we further simplify Equation (35) The final step is the transition from the line integrals to the area integrals over the regions bounded by closed contours of the Montgomery potential (M = M 0 ): Equation (37) shows that there could be no net diapycnal transport across broad regions on the isopycnal surfaces laterally bounded by closed streamlines, which confirms and rationalizes the double-diffusive insulation effect discussed and numerically modeled in Section 3. The foregoing arguments are certainly not relevant for the turbulent ocean.In the turbulent case, the flux laws do not satisfy Equation ( 8) and therefore Equation ( 30)-and all that follows-does not apply.Hence, no analogous constraints on diapycnal velocity are expected to arise in purely turbulent systems.

Discussion
This communication draws attention to the dramatic differences in the dynamics of diapycnal transport realized when the dominant mixing agent is double-diffusion and when mixing is controlled by mechanically-generated turbulence.Double-diffusion occurs in numerous geophysical and astrophysical fluid systems.In the ocean, nearly half of water masses are potentially affected by double-diffusive mixing [42].Double-diffusion is also likely to play a major role in the dynamics of massive stars and giant planets [43,44].Unlike mechanically-generated turbulence, double-diffusion is characterized by unequal eddy diffusivities of density components, which depend on the background density ratio.This feature is ultimately responsible for rather unique properties of double-diffusive mixing, such as its ability to actively constrain diapycnal volume transport-all smooth steady solutions of the T-S advection-diffusion equations are necessarily characterized by low diapycnal velocities.This property is contrasted with the passive behavior of systems dominated by mechanically generated turbulence, where it is possible to construct consistent steady-state solutions for an arbitrarily wide range of diapycnal velocities.
Our interest in the problem of diapycnal volume transport was ignited by the need to fully understand dynamics of the oceanic thermohaline circulation and associated meridional overturning.The classical view [9] ascribes small-scale diapycnal mixing a critical role for its establishment and maintenance.This notion has profoundly affected the evolution of physical oceanography.It motivated development of the extensive small-scale mixing research program, largely aimed at quantification of the average diapycnal diffusivity in the ocean [16,45].It is commonly assumed that the increase in diapycnal diffusion of temperature and salinity necessarily amplifies the diabatic component of the meridional overturning [46].Our study presents a peculiar counter-example of this tendency by showing that mixing can have an adverse impact on overturning.Using a suite of basin-scale simulations and theoretical models, we demonstrate that the inclusion of double-diffusion can place rather severe constraints on the magnitude of vertical diapycnal velocity.In the extreme case when vertical mixing is dominated by double-diffusion, typical values of diapycnal velocity (w * ∼ 5 × 10 −9 m/s) are less, by at least an order of magnitude, than w * driven by mechanically-generated turbulence with comparable T-S diffusivities.Simulations indicate that these constraints are realized both locally and in the isopycnal-average sense.In essence, double-diffusion acts to seal the thermocline by preventing the leakage of seawater (both upward and downward) across the isopycnal surfaces on which double-diffusion is a dominant mixing process.Mechanically generated turbulence, on the other hand, offers no restrictions on diapycnal velocity.When both double-diffusion and turbulence are present, the constraints on w * loosen (tighten) with the increase (decrease) of the fraction of the overall mixing attributed to turbulence.
When the contributions to mixing from double-diffusion and turbulence are comparable-the regime which is perhaps realized in the central Atlantic thermocline [32]-the maximum diapycnal velocity is on the order of ∼ 10 −7 m/s.This value is comparable to current estimates of diapycnal velocity, which raises an intriguing possibility that diapycnal transport could be controlled by constraints on w * imposed by double-diffusion.This suggestion is distinct from the ideas expressed by mainstream models of thermohaline circulation [47] in which the T-S advective-diffusive dynamics constitute only a part of the problem.Ultimately, diapycnal velocity is selected by invoking three-dimensional large-scale balances involving both momentum and density.It should be realized, however, that the latter proposition is based on mixing models with uniform and equal vertical diffusivity (K T = K S = const).The examples in this paper indicate that such models may not capture all the relevant dynamics and therefore we urge caution in conceptualizing the thermohaline circulation based on oversimplified mixing parameterizations.
Our attempts to force the double-diffusive thermocline by applying a large vertical velocity which lies outside of the permitted range-the range in which smooth regular solutions can be found-have led to rather dramatic transformations of the stratification into a series well-defined convective layers separated by stratified interfaces.In this regard, it should be noted that one of the most interesting and controversial problems of double-diffusive theory concerns the origin and dynamics of so-called thermohaline staircases-stepped structures in the vertical temperature and salinity profiles.Several hypotheses have been proposed to explain the formation mechanisms [39,48] and our analysis of the "forced" experiments (Section 3.5) may add another candidate to the list.Thermohaline staircases could form as a response to strong externally imposed vertical advection, provided that the magnitude of vertical velocity is such that the advective-diffusive balance of Munk's [9] type can only be satisfied by discontinuous step-like solutions.
Of course, there are a number of uncertainties in the presented model, particularly with regard to its ability to incorporate all relevant dynamics of the oceanic circulation.For instance, it is not clear how resilient the double-diffusive insulating blanket is in the presence of active mesoscale variability, which can also impact diapycnal transport [38].The processes discussed in this study could be affected by the nonlinearities of the equation of state [49].Thus, the quantitative accuracy of double-diffusive insulation theory is readily questioned.However, the differences in the way double-diffusion and turbulence affect diapycnal transport identified here are suggestive and are likely to be realized in Nature.In this regard, it should be emphasized that large-scale modeling studies which take into account double-diffusion [26,27,50] report systematic reduction in the strength of meridional overturning.Although the extent of this reduction varies considerably between models, it is possible that these results can be attributed to the double-diffusive insulation effect.Finally, it should be mentioned that while the specific solutions in this paper are based on the assumed parameterizations of double-diffusive fluxes, which require further refinements and testing, additional experiments (not shown) indicate that our results are not particularly sensitive to the choice of the functional relations for K T (R ρ ) and K S (R ρ ).Therefore, we expect our qualitative conclusions to be sufficiently robust.
where ε << 1.As is common in weakly nonlinear models [51], we search for a solution of the governing Equation (39) using the power series in ε: T = T 0 (z) + εT 1 (z) + ε 2 T 2 (z) + ... S = S 0 (z) + εS 1 (z) + ε 2 S 2 (z) + ... w = w 0 + εw 1 + ε 2 w 2 + ... (41) These power series are substituted in the governing Equation (39) and the terms of the same order in ε are collected.At the zero order, we essentially arrive at the constant flux ratio model discussed in Section 2.2, which requires zero diapycnal velocity.The zero order equations for T 0 (z) and S 0 (z) are satisfied by the linear gradients: The boundary conditions for the individual T-S components become: (1) for i = 1, 2, ... .
The first order equations are given by Multiplying the second equation by γ 0 and subtracting from the first yields: and integrating Equation (44) twice in z, subject to Equation (45) and the boundary conditions (Equation ( 43)), we arrive at: Using Equations ( 45) and ( 46), we write down the second order equations: Multiplying the second equation in Equation ( 47) by γ 0 and subtracting from the first one yields: The third order T-S equations are obtained in a similar manner.When the third order salinity equation is multiplied by γ 0 and subtracted from the third order temperature equation, we obtain: and third order respectively.The analytical curves rapidly converge to the numerical solution with increasing order of the expansion.The third order asymptotic closely approximates the fully nonlinear solution even for the relatively large ε used for the calculation in Figure A1.
the range of diapycnal velocities realized in three-dimensional double-diffusive simulations (Section 3) is of the same order as Equation ( 60) and much less than in the corresponding turbulent experiments (e.g., Figures 7 and 8).T is computed numerically (black curve) and compared with the calculation based on the asymptotic expansion in ε truncated at the first, second, and third orders (blue, green, and red curves respectively).For ε = 0.32 (δ = 0.1), the numerical result is almost indistinguishable from the third order asymptotic prediction.

Fluids 2016, 1 , 25 4 of 27 Figure 1 .
Figure 1.One-dimensional model.We search for the temperature and salinity profiles ( ) T z and ( ) S z satisfying the vertical advection-diffusion equations for given vertical velocity ( w ) and boundary conditions at the ends of the mixing zone (

.
It should be noted that the advective-diffusive balance (Equation (2)) remains passive even when the diffusivities of heat and salt are spatially uniform but not equal to each other (

Figure 1 .
Figure 1.One-dimensional model.We search for the temperature and salinity profiles T(z) and S(z) satisfying the vertical advection-diffusion equations for given vertical velocity (w) and boundary conditions at the ends of the mixing zone (−H bot < z < −H top ).

Fluids
a quantity measuring the relative contributions from turbulence and double-diffusion to the net mixing.

Figure 2 .R
Figure 2. The range of diapycnal velocities permitted in the hybrid model which includes both double-diffusive and turbulent mixing.Parameter ( ) 0

Figure 2 .
Figure 2. The range of diapycnal velocities permitted in the hybrid model which includes both double-diffusive and turbulent mixing.Parameter K turb /K dd T (R 0 ) measures the relative contributions from turbulence and double-diffusion to the net mixing.The numerical calculations resulting in regular solutions are indicated by heavy dots and light dots represent conditions under which no solutions were found.

Fluids
The model solves the incompressible Boussinesq equations of motion with the linear equation of state.The model geometry is a rectangular domain with dimensions (L x , L y , L z ) = (6400 km, 6400 km, 3000 m) (16) resolved by 64 × 64 × 60 Cartesian grid points with ∆x, ∆y = 10 5 m and ∆z = 50 m.The relatively coarse horizontal resolution prevents generation of mesoscale variability, thereby isolating phenomena of interest-double-diffusive and turbulent diapycnal mixing.Thermohaline forcing is represented by strong relaxation of the surface temperature and salinity to the zonally uniform target patterns shown in Figure 3a,b.The surface temperature and salinity vary linearly with latitude from T South = 30 • C and S South = 37.75 • C to T North = 5.67 • C, and S South = 34.22• C. Surface forcing also includes y-dependent zonal wind stress (Figure 3c) prescribed as follows:

Fluids 2016, 1 Figure 3 .
Figure 3.The meridional patterns of the model forcing fields.Thermohaline forcing is applied by relaxing the surface temperature and salinity to the target patterns shown in (a) and (b) respectively.The wind stress is shown in (c).All forcing fields are zonally-uniform.
to maintain the numerical stability of the model.

Figure 3 .
Figure 3.The meridional patterns of the model forcing fields.Thermohaline forcing is applied by relaxing the surface temperature and salinity to the target patterns shown in (a) and (b) respectively.The wind stress is shown in (c).All forcing fields are zonally-uniform.

Figure 4 .
Figure 4.The final state ( 200 years t = ) realized in the numerical experiments with double-diffusive mixing (the constant flux ratio case): (a) The horizontal temperature distribution and velocity pattern at 500 m z = − ; (b) The zonal section of temperature at 0.5 y y L = .Clearly visible is a well-defined thermocline with relatively warm water extending several hundred meters downward from the surface.

Figure 5 .
Figure 5. Three-dimensional view of the average isopycnal surface represent the slopes of isopycnal surfaces in x and y.For steady-state circulation patterns, Equation (20) reduces to

Figure 4 .
Figure 4.The final state (t = 200 years) realized in the numerical experiments with double-diffusive mixing (the constant flux ratio case): (a) The horizontal temperature distribution and velocity pattern at z = −500 m; (b) The zonal section of temperature at y = 0.5L y .Clearly visible is a well-defined thermocline with relatively warm water extending several hundred meters downward from the surface.

Figure 4 .
Figure 4.The final state ( 200 years t = ) realized in the numerical experiments with double-diffusive mixing (the constant flux ratio case): (a) The horizontal temperature distribution and velocity pattern at 500 m z = − ; (b) The zonal section of temperature at 0.5 y y L = .Clearly visible is a well-defined thermocline with relatively warm water extending several hundred meters downward from the surface.

Figure 5 .
Figure 5. Three-dimensional view of the average isopycnal surface represent the slopes of isopycnal surfaces in x and y.For steady-state circulation patterns, Equation (20) reduces to

Figure 5 .
Figure 5. Three-dimensional view of the average isopycnal surface ρ = ρ av defined in Equation (19).Color coding represents the density ratio distribution on this isopycnal.The observed range of density ratios 1 < R ρ < 3 indicates that fingering is active and that its intensity is relatively high.

Figure 6 .
Figure 6.Physical interpretation of diapycnal velocity.A Lagrangian particle initially located on a motionless tilted isopycnal surface is advected from A to B by a cross-isopycnal flow.Point C is the vertical projection of B onto the isopycnal surface.Diapycnal velocity represents the rate of the increase in vertical separation of the Lagrangian point from the isopycnal surface.

Figure 6 .
Figure 6.Physical interpretation of diapycnal velocity.A Lagrangian particle initially located on a motionless tilted isopycnal surface is advected from A to B by a cross-isopycnal flow.Point C is the vertical projection of B onto the isopycnal surface.Diapycnal velocity represents the rate of the increase in vertical separation of the Lagrangian point from the isopycnal surface.
,b becomes particularly striking when we recall that typical T-S diffusivities in the two experiments are comparable.For instance, the volume-averaged diffusivities in the double-diffusive experiment are K dd T , K dd S = (1.07, 2.18) × 10 −5 m 2 /s, and these values are similar to the turbulent diffusivity (Equation (26)).

Figure 8 ..
Figure 8.The vertical distribution of diapycnal transport.(a) The mean diapycnal velocity evaluated at various isopycnals ( * is o w ) and plotted as a function of the average depth of those surfaces; (b) The local vertical profile of diapycnal velocity ( * lo c w ) at ( ) ( ) , 0.5 ,0.5 x y x y L L = .Both diagnostics indicate that diapycnal velocity in the turbulent case (indicated by the blue curves) substantially exceeds that in the double-diffusive ocean (green curves).The patterns of * i s o w and * lo c w are qualitatively

Figure 9 K 9 ~5
presents the variation of the isopycnal-averaged diapycnal velocity * 19) as a function of turb K .These diagnostics also support our interpretation of double-diffusive mixing as an insulating mechanism.When turb K , the overall system behaves like the fully double-diffusive case (e.g., Figure7a).Values of mean diapycnal velocity are on the order of

Figure 8 .
Figure 8.The vertical distribution of diapycnal transport.(a) The mean diapycnal velocity evaluated at various isopycnals (w * iso ) and plotted as a function of the average depth of those surfaces; (b) The local vertical profile of diapycnal velocity (w * loc ) at (x, y) = 0.5L x , 0.5L y .Both diagnostics indicate that diapycnal velocity in the turbulent case (indicated by the blue curves) substantially exceeds that in the double-diffusive ocean (green curves).The patterns of w * iso and w * loc are qualitatively similar.

Fluids 2016, 1 , 25 13 of 27 Figure 9 .
Figure 9.The variation of the mean diapycnal velocity at the density surface

Figure 9 .
Figure 9.The variation of the mean diapycnal velocity at the density surface ρ = ρ av (w * av ) as a function of K turb .The diagnostics are based on a series of simulations which incorporate both double-diffusive and turbulent mixing.Note the monotonic increase in w * av with increasing K turb .

FluidsFigure 10 .
Figure 10.The assumed patterns of the flux ratio (a) and salt diffusivity (b) used for parameterization of double-diffusion in the numerical simulations (Section 3.4).

8 5
diapycnal velocity calculated for each simulation along the density surface av ρ ρ = -and these results are summarized in Figure 11.The diapycnal velocity for the δ 0 = case is very weak ( * 9 10 m / s w −  ) and negative.However, as δ increases and the flux ratio becomes less and less uniform, * av w monotonically increases.Finally, for δ 1 = , which corresponds to the turbulent flux ratio model ( ) γ R ρ =, the average diapycnal velocity increases to * case become comparable to those realized in the fully turbulent model (Figure7b).The pattern of * ( δ ) a v w in Figure11can be reasonably well described as a linear relation and is

Figure 10 .
Figure 10.The assumed patterns of the flux ratio (a) and salt diffusivity (b) used for parameterization of double-diffusion in the numerical simulations (Section 3.4).
av w -the isopycnal-averaged diapycnal velocity calculated for each simulation along the density surface av ρ ρ = -and these results are summarized in Figure 11.The diapycnal velocity for the δ 0 negative.However, as δ increases and the flux ratio becomes less and less uniform, * av w monotonically increases.Finally, for δ 1 = , which corresponds to the turbulent flux ratio model ( )

Figure 11 .
Figure 11.A series of large-scale simulations in which the variation in the flux ratio (as measured by the parameter δ) is systematically increased.For each experiment, the mean diapycnal velocity ( * av w ) at

Figure 12 .
Figure 12.Vertical T-S profiles obtained with the "forced" model, in which upward vertical velocity was prescribed below the surface and above the sea-floor of the model ocean.The solid (dashed) curves represent the salinity (temperature) profiles: (a) The experiment in which vertical mixing is represented by turbulent diffusion with uniform and equal T-S diffusivities; (b) The analogous experiment performed with double-diffusive mixing.Note the dramatic differences in the resulting stratification patterns.

Figure 12 .
Figure 12.Vertical T-S profiles obtained with the "forced" model, in which upward vertical velocity was prescribed below the surface and above the sea-floor of the model ocean.The solid (dashed) curves represent the salinity (temperature) profiles: (a) The experiment in which vertical mixing is represented by turbulent diffusion with uniform and equal T-S diffusivities; (b) The analogous experiment performed with double-diffusive mixing.Note the dramatic differences in the resulting stratification patterns.
the subtropical ocean on an arbitrarily chosen isopycnal surface 0 ρ ρ = .On this surface, we select a closed temperature contour ( 0 T T = ).Since temperature, salinity and density are assumed to be uniquely related through the linear equation of state, it follows that salinity is also constant along the chosen contours ( 0 S S = ).These contours closely follow the corresponding streamlines of the ideal basic state on which ( ) are conserved exactly.

Figure 13 .
Figure 13.Schematic diagram illustrating the analytical model of double-diffusive insulation.The model suggests that the average diapycnal velocity in the regions bounded by closed streamlines on the density surfaces in the ocean interior (indicated by grey shading) is zero if vertical mixing is double-diffusive but can be finite in the presence of mechanically generated turbulence.

Figure 13 .
Figure 13.Schematic diagram illustrating the analytical model of double-diffusive insulation.The model suggests that the average diapycnal velocity in the regions bounded by closed streamlines on the density surfaces in the ocean interior (indicated by grey shading) is zero if vertical mixing is double-diffusive but can be finite in the presence of mechanically generated turbulence.

Fluids 2016, 1 , 25 24 of 27 Figure A1 .
Figure A1.The non-dimensional departure of temperature from the linear gradient ( ) ( ) T z T z z ′ = − .T ′ is computed numerically (black curve) and compared with the calculation based on the asymptotic expansion in ε truncated at the first, second, and third orders (blue, green, and red curves respectively).For ε 0.32 = ( δ 0.1 =), the numerical result is almost indistinguishable from the third order asymptotic prediction.

Figure A2 .
Figure A2.The range of diapycnal velocities permitted in the one-dimensional model with variable flux ratio ( γ ).Parameter δ measures the extent of variation in γ , with positive (negative) values corresponding to the increasing (decreasing) γ( ) R ρ relation.The numerical calculations resulting in regular solutions are indicated by heavy dots and light dots represent conditions under which no solutions were found.

Figure A1 .
Figure A1.The non-dimensional departure of temperature from the linear gradient T (z) = T(z) − z.T is computed numerically (black curve) and compared with the calculation based on the asymptotic expansion in ε truncated at the first, second, and third orders (blue, green, and red curves respectively).For ε = 0.32 (δ = 0.1), the numerical result is almost indistinguishable from the third order asymptotic prediction.

Fluids 2016, 1 , 25 24 of 27 Figure A1 .
Figure A1.The non-dimensional departure of temperature from the linear gradient ( ) ( ) T z T z z ′ = − .T ′ is computed numerically (black curve) and compared with the calculation based on the asymptotic expansion in ε truncated at the first, second, and third orders (blue, green, and red curves respectively).For ε 0.32 = ( δ 0.1 =), the numerical result is almost indistinguishable from the third order asymptotic prediction.

Figure A2 .
Figure A2.The range of diapycnal velocities permitted in the one-dimensional model with variable flux ratio ( γ ).Parameter δ measures the extent of variation in γ , with positive (negative) values corresponding to the increasing (decreasing) γ( ) R ρ relation.The numerical calculations resulting in regular solutions are indicated by heavy dots and light dots represent conditions under which no solutions were found.

Figure A2 .
Figure A2.The range of diapycnal velocities permitted in the one-dimensional model with variable flux ratio (γ).Parameter δ measures the extent of variation in γ, with positive (negative) values corresponding to the increasing (decreasing) γ(R ρ )relation.The numerical calculations resulting in regular solutions are indicated by heavy dots and light dots represent conditions under which no solutions were found.