Interactive comment on “ A possibility of large scale intrusions generation in the Arctic Ocean under stable-stable stratification : an analytical consideration ” by Natalia Kuzmina

The author addresses an issue of destabilization of the geostrophic flow by vertical diffusion of density and speculates that such process might be the relevant in explaining the formation of intrusive layers in the Arctic Ocean. I enjoyed reading this paper and it is particularly pleasing to see the theoretical analysis. While I believe that there is a potentially great scientific merit to this work, I have several concerns about the clarity of the presentation, novelty of the theoretical results, and applicability to the Arctic Ocean. Please find my detailed comments and questions below that are aimed at improving the quality of the paper.


Introduction
Study of the intrusions in oceanic frontal zones is necessary to analyse the mechanisms of ventilation and mixing in the ocean interior (see, for example, Zhurbas et al., 1983Zhurbas et al., , 1987Rudels et al., 1999Rudels et al., , 2009Kuzmina and Zhurbas, 2000;Walsh and 20 Ruddick, 2000;Merryfield, 2000;Radko, 2003;Richards and Edwards, 2003;Kuzmina et al., , 2011Smyth 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 double diffusion (Stern, 1967;Ruddick and Turner, 1979;Toole and Georgi, 1981;McDougall, 1985aMcDougall, , 1985bNiino, 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 25 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 horizontal scale reaching more than 100 km (Rudels et al., 1999(Rudels et al., , 2009Kuzmina et al., 2011) observed at stable-stable stratification (i.e., when the mean salinity increases with depth while the mean temperature decreases with depth). It can be suggested that 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 3D interleaving 30 model, taking into account differential mixing at a no baroclinicity front and observations of intrusive layering at a pure thermohaline front in the PDW. Merryfield's (2002) findings were confirmed by Kuzmina et al. (2014). However, the 2D model of interleaving driven by differential mixing at the baroclinic front did not show a satisfactory fit simultaneously between the three modelled parameters, namely the vertical scale, growth time and slope of the fastest growing mode, and observations of intrusions in a frontal zone with substantial baroclinicity in the upper PDW layer (Kuzmina et al., 2014). In particular, it was 35 found that the vertical scale of the most unstable mode is two to three times smaller than the vertical scale of intrusions observed Ocean Sci. Discuss., doi:10.5194/os-2016Discuss., doi:10.5194/os- -15, 2016 Manuscript under review for journal Ocean Sci. Published: 10 May 2016 c Author(s) 2016. CC-BY 3.0 License.
2 in the baroclinic front. Furthermore, it is worth noting that the 2D models of double-diffusive interleaving, as applied to typical baroclinic fronts in the ocean, are able to forecast intrusive layers of no more than 10 m vertical length scale (Kuzmina and Rodionov, 1992;Kelley, 1997, 2001;Kuzmina and Zhurbas, 2000;Kuzmina and Lee, 2005;. Therefore, despite the fact that there are proven-by-simulation hypotheses of a merger of small vertical-scale intrusions into larger structures (Radko, 2007), new approaches to the mathematical description of the formation of large intrusions in areas of 5 baroclinic fronts appear relevant.
We suggest that interleaving at baroclinic fronts may be considered as a result of 3D 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 was studied theoretically by Miles (Miles, 1965). Based on an analogy between the equations describing the dynamics of large-scale atmospheric 10 perturbations and the Orr-Sommerfeld equation (Lin, 1955), Miles analysed the instability of the critical layer (the very thin layer in which the phase velocity of disturbance is equal to the velocity of zonal flow). As a result, Miles built an analytical asymptotic solution taking into account the very small, but finite vertical diffusion of buoyancy. Based on the analysis, Miles concluded that the influence of vertical diffusion of buoyancy in destabilising the zonal wind is not essential in comparison with baroclinic instability (the generation of cyclones and anticyclones) for typical atmospheric geostrophic wind. One can assume, 15 however, that other situations can be observed in the deep ocean. Indeed, in the Polar zones, for example, in the Eurasian Basin of the Arctic, very weak geostrophic currents are observed at deep layers (Aagaard, 1981). These currents can have a large horizontal (transverse) scale and large time scale of variability, the latter being estimated at much more than one 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, it seems reasonable to suggest that the role of diffusion of buoyancy in the destabilization of weak geostrophic 20 currents can be important. Therefore, in such circumstances one would expect the formation of intrusions, rather than vortices.
The present work is devoted to the search for analytical unstable (increasing with time) and stable (decreasing with time) solutions based on the potential vorticity equation describing the 3D dynamics of a weak baroclinic front, with the vertical diffusion of buoyancy included. The results, hopefully, will make it possible to obtain some new conceptions about the causes of the formation of large intrusions, particularly in the regions of the Arctic Basin with stable-stable stratification. 25

Problem formulation, derivation of basic equation, and solution search
Let us consider the problem of the 3D instability of the baroclinic front on the basis of the linearized equations of motion in quasi-geostrophic approximation (see, for example, Pedlosky, 1979;Cushman-Roisin, 1994): is a characteristic value of the buoyancy frequency in the frontal zone, and s and s are dimensional constants, either positive or negative, that characterize the cross-front gradients 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 2D (see, for example, Kuzmina, Rodionov, 1992;May, Kelley, 1997;Kuzmina, Zhurbas, 2000) and 3D (Eady, 1949;Miles, 1965;Smyth, 15 2008) instabilities of the oceanic baroclinic fronts. However, the oceanic fronts can be characterized by not only the cross-front gradient of density, but also the cross-front gradient of the buoyancy frequency. This is the case described by Eq. (5) Indeed, if the mean density distribution is expressed by Eq. (5), the geostrophic current velocity will be Equation (7)  Ocean Sci. Discuss., doi:10.5194/os-2016-15, 2016 Manuscript under review for journal Ocean Sci. Published: 10 May 2016 c Author(s) 2016. CC-BY 3.0 License.

4
As 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 long-wave disturbances (perturbation on a planetary scale) of weak geostrophic current ) which satisfy the following relationship between the vertical and horizontal length scales (H and L , respectively): The relative error committed in the solution by so doing is expected to be of the order of 2  , and the smaller 2  , the smaller the error.
In accordance with our approximation, Eq. (8) corresponds to the advection equation (4) at Thus, the vorticity equation (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, to which this theory is applied is approximately equal to 100 20  H m, the ratio of L U/ is estimated as s -1 . Based on the latter estimate, one can suggest that the vertical circulation caused by the temporal change of vorticity and frictional force will not significantly affect the dynamics of large-scale disturbances. This hypothesis will 20 be tested a posteriori by analysing the solutions obtained.
Given that const w  0 , we can put 0 0  w in Eq. (9), and therefore rewrite Eq. (9) as Based on the 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). 25 Let us now pay attention to an important issue. Namely, if we suppose that 0 ) (  z U in Eq. (10) and consider salt fingering instead of diffusion of buoyancy, then, in addition to Eq. (2), it will be necessary to write the following two equations instead of Eq. (10): Ocean Sci. Discuss., doi: 10.5194/os-2016-15, 2016 Manuscript under review for journal Ocean Sci. Equations (10´) along with (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, 5 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. However, the derivation of the model equations, which we have done above, is useful to understand the limits of the model's applicability.
From the point of view of the author of this work, a simple quasi-stationary (geostrophic) system of equations accurately describes the large scale movement especially in the Arctic Ocean, where the influence of the beta effect is not significant, the 10 baroclinic fronts of large width in the ocean interior are often not intense (Kuzmina et al., 2011), and the baroclinic radius of deformation, f HN / , at 100 H m, does not exceed 2-5 km (see also Section 3).
To analyse the instability of the geostrophic flow, let us consider a layer with a vertical scale of 0 2H and place the coordinate system on the middle line of the layer. For the analysis of the instability in the frame of Eqs. (2) and (10), we will consider a symmetric relative to the midline of the layer geostrophic flow with quadratic z-dependence of velocity: 15 A parabolic dependence of the geostrophic flow velocity upon the vertical co-ordinate 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 not unlikely to observe changes of the buoyancy frequency in the cross-front direction, which indicate the existence of linear shear of geostrophic velocity. Consideration of geostrophic flow instability with the velocity profile of is also 20 possible on the basis of analytical methods, and we will consider the related issues below.
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 require the vanishing of vertical velocity at the layer boundaries. In accordance with our approximation, this condition is satisfied.
Due to the fact that the model takes into account vertical diffusion, it is logical to take the conditions of zero buoyancy flux 25 (for density perturbations) at the layer boundaries: 0 (the type 1 boundary conditions). It is reasonable to consider another type of condition too, namely, the slippery boundary conditions or equivalent conditions of the vanishing density disturbances at the boundaries: (the type 2 boundary conditions). Under the type 2 boundary conditions, it is necessary to require the absence of convergence or divergence of buoyancy flux within the layer: . This requirement is necessary because the convergence or divergence of the buoyancy flux within 30 the layer may increase or, conversely, decrease the stability of the layer.

6
To analyse the instability of geostrophic flow, we will seek the solution of Eq. (11) where k is the wave number along the x axis and с is the growth rate. Disturbances of the horizontal velocities will be expressed as . The solution will be unstable, i.e., increasing with time, if the imaginary part of c is positive: Substituting (12) into (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 ? 10 It is easy to verify that the following functions are the partial solutions of Eq. (13): where for any point 0 z z  lying inside the layer. To test partial solutions (14), one has to substitute ) ( 1 z F and ) ( 2 z F from (14) into Eq. (13), reduce the latter to a cubic 15 and evaluate the coefficients 0 A , 1 A , 2 A , and 3 A . It is easy to obtain that this polynomial is identically zero (i.e., Based on the theory of ordinary differential equations (see, for example, Polyanin and Zaitsev, 2001)   Let us now consider the unstable and stable solutions of (15). 25

Unstable solutions
Solution (15) According to (16) . Under such conditions, the requirement of the absence of the buoyancy flux convergence/divergence within the layer is satisfied: in the case of a flow that is symmetric relative to the 20 midline of the layer and the parity of function (15), 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 satisfied: Obviously, in this case, as in the case of (17) . To make it, we have to seek solutions of (13) (preliminarily rewriting this equation for ) in the form: is constructed similarly to previous considerations (see formula (15) With regard to the boundary conditions, the problem is reduced to an eigenvalue problem. Here, we briefly discuss only the 15 evaluation of the growth rate 2 c . It can be easily shown that the solution ) (z . This equality can be satisfied if n is zero or an imaginary number. The imaginary value of n indicates that the value of the growth rate 2 c is negative for all values of the wave number k . In this case, the disturbance is described by trigonometric functions that decrease with time.

Stable solutions 20
Stable solutions of Eq. (13)  . To construct own functions of the eigenvalue problem for the case 1 The solutions describe slow time-decaying, long waves that can move, in contrast to the Rossby waves, not only to the west but also to the east according to the sign of s (see Eq. (7)). Moreover, if Comparing (16) and (19), we can conclude that the phase velocity of the stable and unstable disturbances has a different sign. That is, stable and unstable perturbations described by solutions (15) will move in different directions with respect to the flow and a fixed observer.
For case 1and type 2 boundary conditions, a plot of the stable solutions corresponding to the disturbances of density is presented in Fig. 3. 5

Obtained solutions: some comments
Our own functions, obtained in the previous subsections, have the vertical structure of the unstable perturbations, which differs significantly from those of the classical 3D and 2D interleaving models for double diffusive interleaving at the oceanic front (Stern, 1967;Toole and Georgi, 1981;McDougall, 1985aMcDougall, , 1985bNiino, 1986;Yoshida et al., 1989;Kuzmina and Rodionov, 1992;May and Kelley, 1997;Kuzmina, 2000;Kuzmina and Zhurbas, 2000). Given that the intrusions in the oceanic fronts have 10 different forms, the present results may be useful for interpreting empirical data. However, the simplicity of our model does not allow us 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 double diffusion 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 neglect of vertical friction. A similar feature is typical of 15 a simple model described above. The growth rate increases with the increase in the wave number k up to the limit k for which the constraint of ) 2 /(  .The formula relating the growth rate and the vertical scale of the disturbances will be the same as the previous one, but with H instead of 0 H . 30 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, if we take into account that in our model . This allows the squared parameter  to be presented as 0 4 0 is a diffusion analogue of the Reynolds number or Peclet number. It is worth noting that Eq. (13), being differentiated, corresponds to a simplified form of the Orr-Sommerfeld equation (see e.g. (Lin, 1955;Stern, 1975)) written under the extra-long-wave approximation. However, there are a number of differences between these equations, namely: (a) the destabilising factor in the Orr-Sommerfeld equation is friction, rather than diffusion, 5 (b) the unknown function is a stream function in the vertical plane, and, finally, (c) to analyse instability in the frame of the Orr-Sommerfeld equation it is suggested that viscosity is vanishingly small but finite. For this reason, the disturbances out of the crirical layer are described by the equation of an ideal fluid, but in the region of a thin critical layer the equation of the so-called "viscous regime" is used (see, for example, Miles, 1965).
To obtain the solutions for the "viscous regime", the Orr-Sommerfeld equation is greatly simplified: only the terms 10 describing derivatives of the unknown function of the 4th and 2nd order are considered (Iordanskiy and Kulikovskiy, 1966).
Moreover, the velocity of the flow , where 0 z is a point lying on the midline of the critical layer. In this regard, the solutions obtained in the present work for a parabolic type of geostrophic flow are different from the solutions of the "viscous regime" of the Orr-Sommerfeld equation, which are expressed by Hankel functions of order 1/3 (Lin, 1955). The unstable modes described by solutions (15) cannot be attributed to the instability of the critical layer: see 15 formulas (16) describing the phase velocity of unstable modes. Thus, there are significant differences between the approaches used to study instability in the frame of the Orr-Sommerfeld equation [see also the paper by Miles (Miles, 1965)] and the approach proposed in the present work.
To conclude this section, we note the following.
The instability of the weak geostrophic flow in the frame of the solutions (15) is an oscillatory instability (the growth rate has 20 real and imaginary components). Generally, using interleaving models (Stern, 1967;Tool and Georgy, 1981;McDougal, 1985aMcDougal, , 1985bNiino, 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: 0 Re  c .). The exceptions to this rule are the interleaving models describing the interleaving in the equatorial fronts. In accordance with the modelling efforts (Richards, 1991;Edwards and Richards, 1999;Kuzmina et al., 25 2004;Kuzmina and Lee, 2005), the instability of the equatorial fronts in the scale of intrusive layering is regarded as an oscillatory instability.
The general solution (15) is one of the classes of solutions of Eq. (13). Thus, for example, at it is possible to find an analytically general solution too. This solution will have a more complex structure than (15). Detailed analytical consideration of unstable modes based on the analysis of different classes of solutions of Eq. (13) taking into account friction 30 may be a subject for further research. To clearly define the range of applicability of our model, it is worth solving the eigenvalue problem for Eq. (7) for small values of parameter 2  by means of numerical methods. This problem may be the subject of further research too. The analytical solutions found can be used to validate numerical solutions of the eigenvalue problems.
Moreover, the analytical solutions obtained give analytical formulas for own functions, phase velocities and growth/decay rates of disturbances that cannot, as a rule, be found exactly from numerical solutions . 35 Ocean Sci. Discuss., doi:10.5194/os-2016-15, 2016 Manuscript under review for journal Ocean Sci. Published: 10 May 2016 c Author(s) 2016. CC-BY 3.0 License.

Application to thermohaline intrusions in the Eurasian Basin of the Arctic Ocean
It is worth evaluating the time of formation of the large-scale intrusions based on the results of the presented model.
According to Kuzmina et al. (2011) (Kuzmina et al., 2011)], the maximal linear shear can be estimated as  s (1-2)10 -7 m -1 s -1 . Vertical diffusivity K can be estimated in the range of  K 10 -6 -10 3·10 -6 m 2 s -1 (Merryfield, 2002) [see also the following paper where the evaluations of coefficients of diffusivity in the Arctic thermocline were considered (Walsh and Carmack, 2003)]. The typical vertical scale of intrusive layering in the fronts of PDW is approximately 30-40 m (Merryfield, 2002;Kuzmina et al., 2014) In closing 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 the cases considered. According to Eqs. (2) and (10) Ocean Sci. Discuss., doi: 10.5194/os-2016-15, 2016 Manuscript under review for journal Ocean Sci. Published: 10 May 2016 c Author(s) 2016. CC-BY 3.0 License.
In this paper, we investigated analytically the 3D instability of a baroclinic front in quasi-geostrophic, long-wave approximation taking into account the vertical diffusion of buoyancy. It is shown for the first time that not only double diffusion, but the diffusion of buoyancy can cause destabilization of the geostrophic flow. Such instability has to be distinguished from the 2D McIntyre instability (McIntyre, 1970), a type of instability due to flow-dependent fluctuations in turbulent diffusivities 5 (Smyth and Ruddick, 2010) and the 2D baroclinic instability due to double diffusion (Kuzmina and Rodionov, 1992;May and Kelley, 1997;Kuzmina and Zhurbas, 2000;Kuzmina and Lee, 2005).
In contrast to the work by Miles (Miles, 1965), in which it was shown that the influence of vertical diffusion of buoyancy is not essential in comparison with the influence of vorticity change to destabilize the zonal flow, we considered the opposite case, when vertical diffusion of buoyancy can play an important role as a destabilization factor of very weak geostrophic current with 10 linear shear and large cross-frontal scale.
The model we developed can be considered as a modification of Stern's (Stern, 1967). However, instead of analysing 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 analysed taking into account the vertical diffusion of density. This model can be useful for describing stable and unstable disturbances of a planetary scale in the polar regions of the ocean under the stable-stable stratification, particularly in the deep 15 water of the Arctic Basin, where weak baroclinic fronts with a large horizontal (cross-frontal) scale and typical temporal variability period of the order of several years or more are observed, and the beta-effect is negligible.
The stable (decaying with time) solutions are shown to describe long-wave disturbances that, 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 mean flow (parabolic dependence of mean velocity upon vertical co-ordinate) has 20 an action upon the dynamics of disturbances and likewise the  -effect.
Unstable (increasing with time) solutions are used to describe the formation of large-scale intrusions in the areas of baroclinic fronts, which are observed in the Arctic Basin in the regions characterized by an absolutely stable stratification, for example, in the upper layer of the PDW in the Eurasian Basin.
The proposed description of intrusions generation in baroclinic fronts can be considered as a possible alternative mechanism 25 relative to the differential mixing. However, at the moment this is just a hypothesis, and further efforts, both in theoretical modelling and field data analysis, are needed to justify it.