Generation of large-scale intrusions at baroclinic fronts: an analytical consideration with a reference to the Arctic Ocean

Analytical solutions are found for the problem of instability of a weak geostrophic flow with linear velocity shear accounting for vertical diffusion of buoyancy. The analysis is based on the potential-vorticity equation in a long-wave approximation when the horizontal scale of disturbances is considered much larger than the local baroclinic Rossby radius. It is hypothesized that the solutions found can be applied to describe stable and unstable disturbances of the planetary scale with respect, in particular, to the Arctic Ocean, where weak baroclinic fronts with typical temporal variability periods on the order of several years or more have been observed and the β effect is negligible. Stable (decaying with time) solutions describe disturbances that, in contrast to the Rossby waves, can propagate to both the west and east, depending on the sign of the linear shear of geostrophic velocity. The unstable (growing with time) solutions are applied to explain the formation of large-scale intrusions at baroclinic fronts under the stable–stable thermohaline stratification observed in the upper layer of the Polar Deep Water in the Eurasian Basin. The suggested mechanism of formation of intrusions can be considered a possible alternative to the mechanism of interleaving at the baroclinic fronts due to the differential mixing.


Introduction
The study of intrusions in oceanic frontal zones is required to understand the mechanism of ventilation and mixing in the ocean interior (see, e.g., Zhurbas et al., 1993Zhurbas et al., , 1987;;Rudels et al., 1999Rudels et al., , 2009;;Kuzmina and Zhurbas, 2000;Walsh and Ruddick, 2000;Merryfield, 2000;Radko, 2003;Richards and Edwards, 2003;Kuzmina et al., 2005Kuzmina et al., , 2011;;Smyth and Ruddick, 2010).Intrusive layering, as a rule, results from the instability of oceanic fronts.One of the major mechanisms responsible for the instability of both thermohaline and baroclinic fronts is related to the double diffusion (Stern, 1967;Ruddick and Turner, 1979;Toole and Georgi, 1981;McDougall, 1985a, b;Niino, 1986;Yoshida et al., 1989;Richards, 1991;Kuzmina and Rodionov, 1992;May and Kelley, 1997;Kuzmina, 2000).However, in the Eurasian Basin of the Arctic Ocean there are baroclinic and thermohaline fronts within the upper layer of the Polar Deep Water (PDW) populated with intrusive layers of vertical length scale as large as 30 m and with horizontal scale reaching up to more than 100 km (Rudels et al., 1999(Rudels et al., , 2009;;Kuzmina et al., 2011) observed at the stable-stable stratification (i.e., increasing for the mean salinity while decreasing for the mean temperature with depth).It can be suggested that the thermohaline intrusions within the upper layer of PDW are driven by differential mixing.Merryfield (2002) was the first to show satisfactory agreement between calculations of unstable modes from a three-dimensional (3-D) interleaving model that accounted for the differential mixing at a non-baroclinic front and observations of intrusive layering at a pure thermohaline front in the PDW.Findings by Merryfield (2002) were confirmed by Kuzmina et al. (2014).However, the 2-D model of interleaving driven by differential mixing at the baroclinic front failed to simultaneously fit three modeled parameters, namely, the vertical scale, the growth time, and the slope of the fastest growing mode, with observations of intrusions in a frontal zone with a substantial baroclinicity in the upper PDW layer (Kuzmina et al., 2014).In particular, it was found that the vertical scale of the most unstable mode was about 2 to 3 times smaller than the vertical scale of intrusions observed in the baroclinic front.Furthermore, it is worth not-Published by Copernicus Publications on behalf of the European Geosciences Union.
N. Kuzmina: Generation of large-scale intrusions at baroclinic fronts ing that the 2-D models of double-diffusive interleaving, as applied to typical baroclinic fronts in the ocean, are able to forecast intrusive layers with vertical length scale of no more than 10 m (Kuzmina and Rodionov, 1992;May andKelley, 1997, 2001;Kuzmina and Zhurbas, 2000;Kuzmina and Lee, 2005;Kuzmina et al., 2005).Therefore, despite the proven-by-simulation hypothesis of intrusions of small vertical scale merging into larger structures (Radko, 2007), new approaches to the mathematical description of the formation of large intrusions in the area of baroclinic fronts become relevant.
We suggest that the interleaving at a baroclinic front may be considered as a result of 3-D instability of weak geostrophic current due to the combined effects of vertical shear and diffusion of density (buoyancy).
The effect of vertical diffusion of buoyancy on the baroclinic instability of geostrophic zonal wind has been studied theoretically by Miles (1965).Proceeding from the analogy between the equations describing the dynamics of large-scale atmospheric perturbations and the Orr-Sommerfeld equation (Lin, 1955;Stern, 1975), Miles (1965) analyzed the instability of the critical layer (a very thin layer in which the phase velocity of a disturbance equals the velocity of zonal flow).This resulted in an analytical asymptotic solution accounting for a very small, though finite, vertical diffusion of buoyancy.Based on the analysis, Miles (1965) concluded that the effect of vertical diffusion of buoyancy on the destabilization of zonal wind is negligible in comparison with the baroclinic instability (the generation of cyclones and anticyclones) for typical atmospheric geostrophic winds.One could assume, however, that other conditions can occur in the deep ocean.Indeed, in the polar zones, for example, in the Eurasian Basin of the Arctic, very weak geostrophic currents have been observed in deep layers (Aagaard, 1981).These currents can have a large horizontal (transverse) scale and large timescale of variability, the latter being estimated to exceed 1 year (Aagaard, 1981).Taking into account that the influence of the β effect on the dynamics of large-scale disturbances is negligible in the Polar Ocean for near the pole, it seems reasonable to suggest that the contribution of diffusion of buoyancy to the destabilization of weak geostrophic currents can be important.Therefore, in such circumstances one would expect the formation of intrusions, rather than vortices.
The present work is devoted to seeking analytical unstable (increasing with time) and stable (decreasing with time) solutions based on the potential-vorticity equation describing the 3-D dynamics of a weak baroclinic front, with the vertical diffusion of buoyancy included.Hopefully the results will provide an opportunity to obtain some new insight into the causes of the formation of large intrusions, particularly in the regions of the Arctic Ocean with the stable-stable stratification.

Problem formulation, derivation of basic equation, and solution search
Let us consider the problem of the 3-D instability of a baroclinic front based on the linearized equations of motion in quasi-geostrophic approximation (see, e.g., Pedlosky, 1992;Cushman-Roisin, 1994): where U and V are zonal and meridional components of the geostrophic velocity; P and ρ are the mean pressure and density both divided by the reference density; N, f , and g are the buoyancy frequency, Coriolis parameter, and gravity acceleration; u, v, and w are velocity fluctuations along the x, y, and z axes, respectively; p and ρ are the pressure and density fluctuations both divided by the reference density; β = ∂f/∂y; = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 ; and the x, y, and z axes are directed eastward, northward and upward, respectively.The vertical friction with a constant coefficient K is considered in the vorticity equation Eq. (3).The density balance equation (Eq.4) takes into account, apart from the advection terms, only the vertical diffusion with a constant coefficient K.The constant coefficients K and K are treated as the average values over the ocean layer under investigation.
Let us take the distribution of mean density, divided by the reference density, as follows: where N 0 = const > 0 is a characteristic value of the buoyancy frequency in the frontal zone, and s and s are dimensional constants, either positive or negative, which characterize the cross-front gradient of density and the vertical shear of the basic geostrophic flow.
The first term on the right of Eq. ( 5) has not been taken into account in the interleaving models describing both the 2-D (see, e.g., Kuzmina and Rodionov, 1992;May and Kelley, 1997;Kuzmina andZhurbas, 2000) and3-D (Eady, 1949;Miles, 1965;Smyth, 2008) instabilities of the oceanic baroclinic fronts.Meanwhile, the oceanic fronts can be characterized by not only the cross-front gradient of density, but also by the cross-front gradient of the buoyancy frequency.This is the case described by Eq. ( 5); the squared buoyancy frequency N 2 = −gdρ/dz is a linear function of y.This dependence is assumed to be weak; |s| f L N 2 0 , where L is the characteristic lateral length scale (width) of the frontal zone (0 ≤ y ≤ L).However, even a weak lateral change in the buoyancy frequency indicates the existence of a quadratic dependence of geostrophic velocity on the vertical co-ordinate z.Indeed, if the mean density distribution is expressed by Eq. ( 5), the geostrophic current velocity will be where U 1 and U 2 are the constituents of geostrophic velocity with linear (U 1 ) and constant (U 2 ) vertical shear: dU 1 /dz = sz, dU 2 /dz = s and U 3 is the barotropic (constant) velocity addition.
As it can be seen from Eq. ( 7), the last term on the left can strengthen or weaken, depending on the sign of s, the impact of the β effect on the dynamics of disturbances.
We will consider at K ≈ K the long-wave disturbances (i.e., perturbations of the planetary scale) of weak geostrophic current that satisfy the following relationship between the vertical and horizontal length scales (H and L, respectively): L L R , where L R = N 0 H /f is the baroclinic Rossby radius of deformation.If we apply Eqs. ( 1)-( 4) to describe the motion in the Arctic Basin, the β-effect term can be ignored because β ≈ 0 in the vicinity of the North Pole.
Taking into account the abovementioned conditions, we may use the method of series expansion at small parameter where Bu is the Burger number (see, e.g., Cushman-Roisin, 1994).For ), it is reasonable to consider only the first term of the series.In this case, we can rewrite the potential-vorticity equation in the simplified form: The introduced-by-the-procedure relative error of the solution is expected to be on the order of δ 2 , and the smaller the δ 2 the smaller the error.According to our approximation, Eq. ( 8) corresponds to the density balance equation Eq. ( 4) for w = w 0 = const: The correspondence between Eqs. ( 8) and ( 9) can be checked by differentiating Eq. ( 9) with respect to z and taking into account that ∂p/∂z = −gρ.Thus, the vorticity equation Eq. ( 3) drops out of consideration.Indeed, given that the diffusivity of mass, K, in the oceanic interior (particularly in the deep water of the Arctic Ocean) probably does not exceed the value of 1×10 −5 m 2 s −1 and the vertical length scale of the intrusions, H , which this theory is applied to, is approximately equal to H ≈ 20-100 m, the ratio of U/ L ∼ K/H 2 is estimated as U/ L < 3 × 10 −8 s −1 .Based on the latter estimate, one can suggest that the vertical circulation caused by the frictional force and temporal change of vorticity will not significantly affect the dynamics of large-scale disturbances.This hypothesis will be tested a posteriori by analyzing the solutions obtained.For further discussion it is important to underline that the geostrophic Richardson number Ri = N 2 0 H 2 / U 2 , where U is the characteristic scale of geostrophic velocity, is much larger than unity for very slow currents.
Given that w 0 = const, we can take w 0 = 0 in Eq. ( 9), and therefore rewrite it as Based on the reasoning above, we can conclude that the slow extra-large-scale disturbances of weak geostrophic flow are described by the quasi-stationary system of Eqs. ( 2) and ( 10).
Let us now pay attention to an important issue.Namely, if we suppose that U (z) = 0 in Eq. ( 10) and consider salt fingering instead of diffusion of buoyancy, in addition to Eq. ( 2), it will be necessary to write the following two equations instead of Eq. ( 10): where S and S are the salinity disturbance and mean, K S and γ = αF T / βF S < 1 are the vertical diffusivity of salinity and the flux ratio for salt finger convection, F T and F S are the vertical fluxes of temperature and salinity, α and β are the temperature expansion and salinity contraction coefficients, respectively.Equation ( 10) along with Eq. (2) constitute the system of equations that was used by Stern (1967) to obtain the polynomial dependence between the growth rate of unstable perturbations, wave numbers and hydrological parameters (see Eq. 4 of Stern, 1967).Therefore, the proposed model, which consists of Eqs. ( 2) and ( 10), can in a certain sense be regarded as an analogue of the model by Stern (1967) for investigating the interleaving on a large horizontal scale.
From the point of view of the author of this paper, a simple quasi-stationary (geostrophic) system of equations accurately describes the large-scale movement especially in the Arctic Ocean, where the influence of the β effect is not significant, the baroclinic fronts of large width in the ocean interior are often not intense (Kuzmina et al., 2011), and the baroclinic radius of deformation, H N/f , at H ∼ 100 m does not exceed 2-5 km (see also Sect. 3).
To analyze the instability of the geostrophic flow in the frame of Eqs. ( 2) and ( 10), let us take a layer with the vertical scale of 2H 0 , move the z axis origin to the middle of the layer, and consider a symmetric relative to the midline geostrophic flow with quadratic z dependence of velocity: A parabolic z dependence of the geostrophic flow velocity can be observed in the rotary flow of the intra-pycnocline vortices, as well as in many other ocean flows.In any case, as mentioned above, in the oceanic frontal zones it is likely to observe changes of the buoyancy frequency in the crossfront direction indicating the presence of linear shear of geostrophic velocity.The consideration of the instability of geostrophic flow with the velocity profile of U = U 1 + U 2 + U 3 is also possible by analytical methods, but this issue falls out of the scope of the present study.
Let us discuss the conditions on the boundaries of the layer in relation to the ocean.Keeping in mind the Eady problem (Eady, 1949), one has to set the vertical velocity vanishing at the layer boundaries.Our approximation meets this condition.
Due to our model accounting for the vertical diffusion, it appears reasonable to accept the conditions of zero buoyancy flux (for density perturbations) at the layer boundaries: p zz = 0 at z = ±H 0 (the type 1 boundary conditions).It is reasonable to consider another type of condition too, namely, the slippery boundary conditions or equivalent of density disturbances vanishing at boundaries: dv/dz = du/dz = ρ = 0 at z = ±H 0 (the type 2 boundary conditions).Under the type 2 boundary conditions, it is necessary to set up the absence of convergence or divergence of buoyancy flux within the layer: p zz (z = H 0 ) = p zz (z = −H 0 ).This condition is necessary, because the convergence or divergence of the buoyancy flux within the layer may increase or, conversely, decrease the stability of the layer.
Using Eq. ( 2), we rewrite Eq. ( 10) as where To analyze the instability of geostrophic flow, we will seek the solution of Eqs. ( 2) and (11) for the lateral boundary conditions of v| y=0 = v| y=L = 0, accordingly to Eq. ( 7), in the form where k is the wave number along the x axis and c is the growth rate.For the positive imaginary part of c, I m(c) > 0, the solution will be unstable, i.e., increasing with time.
The substitution of Eq. ( 12) into Eq.( 11) yields the following equation: We are interested in finding an answer to the following question: is it possible to make certain judgements about the possibility of instability of geostrophic flow in a finite vertical layer, based on the analytical solutions of Eq. ( 13), at some values of parameter c?
It is easy to verify that the following functions are partial solutions of Eq. ( 13): where mc, and U 1 + U 3 − Rec = 0 at an arbitrary point z = z 0 in the layer domain.
Proceeding from the theory of ordinary differential equations (see, e.g., Polyanin and Zaitsev, 2001), due to the linearly independent functions F 1 (z) and F 2 (z), we can express the general solution of Eq. ( 13) for ikc = 5a •K +ikU 3 in the form where B 1 , B 2 , and B 3 are arbitrary constants.It is important to note two facts.First, the functions F 1 (z) and F 2 (z) are even functions, while F 3 (z) is an odd function.Second, despite the singularity of integrands in Eq. ( 15) at z = 0, the function F 3 (z) is differentiable at this point.(The latter becomes evident from the asymptotic analysis of function F 3 (z) for z → 0.) Let us now consider the unstable and stable solutions Eq. (15).

Unstable solutions
According to the expression for parameter c (see Eq. 14), ikc = ik(c 1 + ic 2 ) = 5a • K + ikU 3 , the solution for Eq. ( 15) can be unstable for both s > 0 and s < 0. The real and imaginary parts of kc for the unstable (i.e., growing with time) solutions are kc 2 = 2.5 • |s| kK for s < 0, According to Eq. (16a, b), the unstable solution is realized for Rea < 0, and hence, for any finite wave number k, the function F 1 (z), and all its derivatives increase infinitely if z → ±∞.On the other hand, the function F 3 (z) and all its derivatives vanish if z → ±∞.(It can be seen from asymptotic analysis of the integrals that define F 3 (z), if z → ±∞.)Therefore, to prove the instability in a finite layer, it is necessary to show that F (z) for Rea < 0 is an eigenfunction of the eigenvalue problem with the boundary conditions of type 1 or 2 introduced above.To construct physically correct solutions we will consider the following two cases.In case 1, where the vertical scale of the layer corresponds to our approximation, 2H 0 ∼ H f L/N .In case 2, where the vertical scale of the layer significantly exceeds the vertical scales of the disturbances for which our approximation holds true, H 2H 0 .To satisfy the boundary conditions of type 1 and 2 in case 1, we have to take B 3 = 0 because F 3 (z) is an odd function.The type 1 boundary conditions are reduced to F zz = 0 for z = ±H 0 .Thus, the following equality should be met: Given that 2B 2 /B 1 can have different values, the instability in the framework of the solution for Eq. ( 15) does exist, because in a wide range of typical ocean values of H 0 , s, and K, there is a wave number k 0 f (2N 0 H 0 ) for which Eq. ( 17) is satisfied.The type 2 boundary conditions are reduced to F z = 0 at z = ±H 0 .Under such conditions, the requirement of the absence of the buoyancy flux convergence/divergence within the layer is met; in case of parity of F (z) for B 3 = 0 and for the flow symmetry relative to the midline, the values of buoyancy flux at the boundaries are of the same magnitude and direction (sign).Under the type 2 boundary conditions the following equality should be met: Obviously, in this case, as in the case of Eq. ( 17), there is a wave number k 0 f/(2N 0 H 0 ) for which Eq. ( 18) is satisfied.
Figure 1 presents graphic images of unstable solutions in the form of density disturbances Reρ = Re(dF /dz) = ρ for different boundary conditions for case 1.When building the solutions, typical values of hydrological parameters in relation to the Arctic Basin were used (see Sects.2.3 and 3).
In the case 2, we have to take B 1 = 0, B 2 = 0, and consider F 3 (z) as the solution of the eigenvalue problem.Indeed, F 3 (z) and all its derivatives sharply decrease if z → ±∞, and, consequently, on the boundaries of the large verticalscale layer, the function F 3 (z) and all its derivatives should be infinitesimally small.Therefore, the boundary conditions of type 1 and 2 are satisfied.Indeed, the characteristic vertical scale of the decrease of F 3 (z) at z → ±∞ can be evaluated as h ∼ (ks/K) −1/4 .If H 0 is n times as large as h, the value of |F 3 (z)| at |z| = H 0 will be e n 2 times (!) as small as the maximum value of |F 3 (z)|.Note that the maximum value of geostrophic velocity U max increases with H 0 .However, U max should satisfy the condition of U max k/f 1, which can be rewritten as n 2 (ksK) 1/2 /f 1.It is easy to see that this condition is met in a wide range of k, s, and K for n = 5-10.
A plot of the disturbances of density ρ = Re(dF 3 /dz) corresponding to unstable solutions is presented in Fig. 2. It is worth noting that the function ρ = dF 3 /dz is differentiable at z = 0 likewise the function F 3 (z).
Thus, in accordance with Eq. ( 11), which was obtained from Eqs. ( 1) to (10) for Bu 1, Ri 1 and Pr = 1, the large-scale disturbances can be unstable.Such instability has to be distinguished from the diffusive instability (McIntyre, 1970;Baker, 1971;Calman, 1977), which occurs when Ri < (Pr +1) 2 /4 Pr and is absent at Pr = 1.One of the important distinctions between these two models of baroclinic front instability is that in the present model the disturbances are allowed to have a non-zero slope in the along-front direction, whereas in the model of diffusive instability by McIntyre (1970) the slope is taken to be zero.Therefore, the McIntyre's model and other models in which the term ∂p/∂x, where x is the along-front co-ordinate, in the equations of motions is ignored (McIntyre, 1970;Baker, 1971;Calman, 1977;Munro et al., 2010) can be referred to as the 2-D models.From the mathematical point of view, the models that take into account the along-front slope of disturbances, are much more complicated.Indeed, the analysis of instability in the 2-D models ultimately reduces to finding the roots of a polynomial, depending upon the wave numbers and growth rate.The models, which take into account the along-front slope of disturbances, are reduced to the differential equations with variable, z-dependent coefficients, and such problems can be solved analytically only in rare cases.
Also, it is worth noting that the instability described by Eq. ( 11) is not the critical layer instability analyzed by Miles (1965) for a geostrophic current with constant vertical shear based on a similarity between the potential-vorticity and the Orr-Sommerfeld equations.Indeed, the phase velocity of the unstable disturbances in our model satisfies the inequality U 1 + U 3 − Rec = 0 for any value of the vertical co-ordinate z.As we can see from Eq. ( 16), the phase velocity of unstable disturbances is directed along the geostrophic current and exceeds the maximum velocity of the current.In such a case, in the author's opinion, the most likely is the conversion of the kinetic energy of the main flow into the kinetic energy of disturbances.

Stable solutions
Stable solutions of Eq. ( 13) are realized for Rea > 0. In this case F 1 (z) and all its derivatives vanish for z → ±∞, but F 3 (z) and all its derivatives increase infinitely for z → ±∞.To construct our own functions of the eigenvalue problem for case 1 (2H 0 ∼ H f L/N ), we have to take B 3 = 0.The solutions describe slow time-decay, long waves that can move, in contrast to the Rossby waves, not only to the west but also to the east depending on the sign of s (see Eq. 7).Moreover, if |s| > βN 2 0 /f 2 (which is quite possible especially in polar regions), the long-wave dynamics in the β-plane approximation is determined by the linear shear of geostrophic flow rather than the β effect.The real and imaginary parts of the growth rate of stable perturbations are According to Eq. ( 19), the condition 16) and ( 19), we can conclude that the phase velocity has a different sign for the stable and unstable disturbances.That is, the stable and unstable perturbations described by solutions for Eq.(15) will move in opposite directions with respect to the flow and a fixed observer.
For case 1 type 2 boundary conditions, a plot of the density disturbances Reρ = Re(dF /dz) = ρ corresponding to the stable solutions is presented in Fig. 3.

Obtained solutions: discussion
Our model does not allow one to determine the maximum growth rate.Here again we can see an analogy with the work by Stern (1967).Indeed, in a well-known paper by Stern (1967), which was the first study of the doublediffusion instability of the infinite thermohaline front, the magnitude of the fastest growing mode was not found.The reason is that the growth rate in Stern's model could indefinitely increase with the horizontal wave number due to the neglected vertical friction.A similar feature is typical for our model.The growth rate increases with the increase in the wave number k up to the limit k for which the constraint of k f/(2H 0 N 0 ) is still valid.Nevertheless, for a rough estimate of the time of formation of unstable perturbations, it is reasonable to use Eq. ( 16).It is also worth evaluating the relationship between the growth rate of unstable disturbances and the layer thickness (case 1) or the characteristic vertical scale of disturbances (case 2).Let us address Eq. ( 17), which follows from the boundary conditions for one of the problems of instability in a finite layer.The parameter χ = Re(−aH 2 0 ) = 0.5(ks/K) 1/2 H 2 0 governs Eq. ( 17).The higher the value of χ , the larger the wave number of the unstable mode for the given values of the problem parameters K, s, and H 0 , and therefore, the larger the growth rate.However, the applicability of our model imposes a constraint on the space of wave numbers, k f/(2N 0 H 0 ).In order to satisfy these two conditions simultaneously in the wide range of variability of hydrological parameters in the ocean, it is reasonable to put 1 ≤ χ ≤ 2. For χ = 2, taking into account Eq. ( 16), we obtain the following formula relating the growth rate of disturbances and the vertical scale of the layer: kc 2 = 1/T = 10 • K/H 2 0 .It is easy to understand the physical meaning of the parameter χ .This parameter characterizes the ratio of advection and vertical diffusion terms depending on the wave number k.Indeed, recalling that in our model U = U 1 + U 3 and the geostrophic velocity on the boundaries of the layer is zero, the maximum velocity at the midline of the layer shall be U max = |s|H 2 0 /2.This allows the squared parameter χ to be presented as χ 2 = 0.25(ks/K)H 4 0 = 0.5 • R d kH 0 , where R d = 0.5sH 3 0 /K = U max H 0 /K is a diffusion analogue of the Reynolds number called the Peclet number.
To conclude this section, we note the following.The instability of the weak geostrophic flow in the frame of the solutions Eq. ( 15) is an oscillatory instability (the growth rate has real and imaginary components).Generally, using interleaving models (Stern, 1967;Toole and Georgy, 1981;McDougal, 1985a, b;Niino, 1986;Yoshida et al., 1989;Kuzmina and Rodionov, 1992;May and Kelley, 1997;Kuzmina and Zhurbas, 2000;Walsh and Ruddick, 2000;Merryfield, 2002), it is possible to obtain the monotonous unstable modes only (the phase velocity of the disturbances is equal to zero: Rec = 0).The exceptions to this rule are the interleaving models related to equatorial fronts.In accordance with the modeling efforts (Richards, 1991;Edwards and Richards, 1999;Kuzmina et al., 2004;Kuzmina and Lee, 2005), the instability of the equatorial fronts in the scale of intrusive layering is regarded as an oscillatory instability.
General solution Eq. ( 15) is one of the classes of solutions of Eq. ( 13).Thus, for example, at ikc = 9a•K+ikU 3 it is also possible to find an analytical general solution of Eq. ( 13).This solution would have a more complex structure than the solution Eq. ( 15).The detailed analytical consideration of unstable modes based on the analysis of different classes of solutions of Eq. ( 13) taking into account the friction may be a subject for further research.In order to clearly define the range of applicability of our model, it would be worth solving the eigenvalue problem for Eq. ( 7) for small values of parameter δ 2 by means of numerical methods.This problem also may be a subject for further research.The analytical solutions found can be used to validate numerical solutions of the eigenvalue problems.Moreover, the analytical solutions obtained provide analytical expressions for eigenfunctions, phase velocities and growth/decay rates of disturbances that cannot, as a rule, be found exactly from numerical solutions.
3 Application to thermohaline intrusions in the Eurasian Basin of the Arctic Ocean It is worth evaluating the time of formation of large-scale intrusions based on the results of the presented model.According to Kuzmina et al. (2011), in the upper layer of the PDW where the large-scale intrusions are observed in the Eurasian Basin at stable-stable stratification, the following estimates of N, f , and β are typical: N ≈ 2 × 10 −3 s −1 , f = 1.4×10 −4 s −1 , and β < 0.3×10 −11 m −1 s −1 (at latitude of 83 • N and higher).Therefore, for disturbances, for example, with the vertical scale of H = 100 m, the Rossby radius of deformation is only H N/f ≈1 km.
According to the derivation of Eq. ( 7), the value of the linear shear s is limited by the inequality of |s| f L N 2 0 .Given that the horizontal scale of the baroclinic fronts (along the cross-front axis y) in the upper layer of the PDW is approximately L ≈ 50-100 km (see examples of transection across the fronts of different types observed in the PDW; Kuzmina et al., 2011), the maximal linear shear can be estimated as |s| ≈ (1 − 2) × 10 −7 m −1 s −1 .Such value of the linear shear is large enough to neglect the β-effect term relative to the linear shear term in Eq. ( 7): βN 2 0 /f 2 s < 0.6 × 10 −2 1.The vertical diffusivity K can be estimated in the range of K = (1−3)×10 −6 m 2 s −1 (Merryfield, 2002;Walsh and Carmack, 2003).We suggest a weak turbulence regime in the layer under consideration: Pr > 1, Pr •Bu 1.The typical vertical scale of intrusive layering in the fronts of PDW is approximately 30-40 m (Merryfield, 2002;Kuzmina et al., 2014).Let us evaluate the time of formation of intrusions with the vertical scale of ∼ 40 m.Using formula k 0 c = 10•K/H 2 0 (see Sect. 2.3), we can estimate the time of formation of the unstable mode as 1/(k 0 c 2 ) for ∼ 5 years for K = 10 −6 m 2 s −1 and ∼ 2 years for K = 3 × 10 −6 m 2 s −1 .
In the closing of this section, let us justify the assumption that the circulations associated with changes in vorticity p are not essential in the description of the formation of intrusions in all considered cases.According to Eqs. ( 2) and ( 10), the characteristic scale of vertical velocity in such circulations can be written as w 1 ∼ U • u • H 0 • k 2 0 /f .In all the considered above cases of application of the model to the Arctic intrusions, the relation of U • k 0 < 10 −8 s −1 is satisfied.Given that small disturbances of horizontal velocity cannot exceed the value of geostrophic velocity U , we find w 1 < 4 × 10 −11 m s −1 .A liquid particle with such vertical velocity travels less than 0.004 m over the period of formation of intrusion (t * = 1/(k 0 c 2 ) ≈ 3 years), while due to the vertical diffusion, the particle displacement is estimated as √ K • t * ≈ 40 m (i.e., 4 orders of magnitude larger).Note also that the decreasing with |z| solution F 3 (z) can be used for the description of generation of intrusions, even if the vertical velocity is not negligibly small.

Conclusions
In this paper, we analytically investigated the instability of a baroclinic front in the quasi-geostrophic, long-wave approximation taking into account the vertical diffusion of buoyancy.Such instability has to be distinguished from the 2-D McIntyre instability (McIntyre, 1970), the instability due to flow-dependent fluctuations in turbulent diffusivities (Smyth and Ruddick, 2010), and the 2-D baroclinic instability due to the double diffusion (Kuzmina and Rodionov, 1992;May and Kelley, 1997;Kuzmina and Zhurbas, 2000;Kuzmina and Lee, 2005).
In contrast to the paper by Miles (1965), who showed that the vertical diffusion of buoyancy is not essential in comparison with the vorticity change in the destabilization of zonal flow, we considered the opposite case, where the vertical diffusion of buoyancy can play an important role as a destabilizer of a very weak geostrophic current with linear shear and large cross-frontal scale.
The model we developed can be considered as a modification of Stern's model (Stern, 1967).However, instead of analyzing the instability of a purely thermohaline front due to the double diffusion (Stern, 1967), in our case the instability of a weak baroclinic front is analyzed taking into account the vertical diffusion of density.This model can be useful for describing stable and unstable disturbances of the planetary scale in the polar regions of the ocean under the stable-stable stratification, particularly in the deep water of the Arctic Ocean, where weak baroclinic fronts with a large horizontal (cross-frontal) scale and typical temporal variability period on the order of several years or more have been observed, and the β effect is negligible.
The stable solutions are shown to describe long-wave disturbances, which, unlike Rossby waves, can move not only to the west but also to the east, depending on the magnitude and sign of the linear shear of geostrophic velocity.It is important to underline that the linear shear of the mean flow (parabolic z dependence of the mean velocity) affects the dynamics of disturbances as well as the β effect.
The unstable solutions can contribute to better understanding of the formation of large-scale intrusions at baroclinic fronts of the Arctic Ocean in the layers characterized by absolutely stable thermohaline stratification, for example, in the upper layer of the PDW in the Eurasian Basin.It is important that the vertical scale of the new modes of instability can reach tens of meters of magnitude, just in accordance with the observations.However, the model is so complex that obtaining the comprehensive results of modeling that can be fully comparable with the empirical data would still remain a future task.