% This is file shear_p1.tex

\NeedsTeXFormat{LaTeX2e}

\documentclass{jfm}
%\documentclass[referee]{jfm} %for double spaced output for submission

% See if the author has AMS Euler fonts installed: If they have, attempt
% to use the 'upmath' package to provide upright math.

\usepackage{graphicx,epsfig}
\usepackage{natbib}
\usepackage{amsmath}
\usepackage{lscape}	%Enables making landscape pages
\usepackage{pdflscape}
\usepackage{color}
\usepackage{lineno} %Used for putting line numbers
\ifCUPmtlplainloaded \else
  \checkfont{eurm10}
  \iffontfound
    \IfFileExists{upmath.sty}
      {\typeout{^^JFound AMS Euler Roman fonts on the system,
                   using the 'upmath' package.^^J}%
       \usepackage{upmath}}
      {\typeout{^^JFound AMS Euler Roman fonts on the system, but you
                   dont seem to have the}%
       \typeout{'upmath' package installed. JFM.cls can take advantage
                 of these fonts,^^Jif you use 'upmath' package.^^J}%
       \providecommand\upi{\pi}%
      }
  \else
    \providecommand\upi{\pi}%
  \fi
\fi

% See if the author has AMS symbol fonts installed: If they have, attempt
% to use the 'amssymb' package to provide the AMS symbol characters.

\ifCUPmtlplainloaded \else
  \checkfont{msam10}
  \iffontfound
    \IfFileExists{amssymb.sty}
      {\typeout{^^JFound AMS Symbol fonts on the system, using the
                'amssymb' package.^^J}%
       \usepackage{amssymb}%
       \let\le=\leqslant  \let\leq=\leqslant
       \let\ge=\geqslant  \let\geq=\geqslant
      }{}
  \fi
\fi

% See if the author has the AMS 'amsbsy' package installed: If they have,
% use it to provide better bold math support (with \boldsymbol).

\ifCUPmtlplainloaded \else
  \IfFileExists{amsbsy.sty}
    {\typeout{^^JFound the 'amsbsy' package on the system, using it.^^J}%
     \usepackage{amsbsy}}
    {\providecommand\boldsymbol[1]{\mbox{\boldmath $##1$}}}
\fi


%%% Example macros (some are not used in this sample file) %%%

% For units of measure
\newcommand\dynpercm{\nobreak\mbox{$\;$dyn\,cm$^{-1}$}}
\newcommand\cmpermin{\nobreak\mbox{$\;$cm\,min$^{-1}$}}

% Various bold symbols
\providecommand\bnabla{\boldsymbol{\nabla}}
\providecommand\bcdot{\boldsymbol{\cdot}}
\newcommand\biS{\boldsymbol{S}}
\newcommand\etb{\boldsymbol{\eta}}

% For multiletter symbols
\newcommand\Real{\mbox{Re}} % cf plain TeX's \Re and Reynolds number
\newcommand\Imag{\mbox{Im}} % cf plain TeX's \Im
\newcommand\Rey{\mbox{\textit{Re}}}  % Reynolds number
\newcommand\Pran{\mbox{\textit{Pr}}} % Prandtl number, cf TeX's \Pr product
\newcommand\Pen{\mbox{\textit{Pe}}}  % Peclet number
\newcommand\Ai{\mbox{Ai}}            % Airy function
\newcommand\Bi{\mbox{Bi}}            % Airy function

% For sans serif characters:
% The following macros are setup in JFM.cls for sans-serif fonts in text
% and math.  If you use these macros in your article, the required fonts
% will be substitued when you article is typeset by the typesetter.
%
% \textsfi, \mathsfi   : sans-serif slanted
% \textsfb, \mathsfb   : sans-serif bold
% \textsfbi, \mathsfbi : sans-serif bold slanted (doesnt exist in CM fonts)
%
% For san-serif roman use \textsf and \mathsf as normal.
%
\newcommand\ssC{\mathsf{C}}    % for sans serif C
\newcommand\sfsP{\mathsfi{P}}  % for sans serif sloping P
\newcommand\slsQ{\mathsfbi{Q}} % for sans serif bold-sloping Q

% Hat position
\newcommand\hatp{\skew3\hat{p}}      % p with hat
\newcommand\hatR{\skew3\hat{R}}      % R with hat
\newcommand\hatRR{\skew3\hat{\hatR}} % R with 2 hats
\newcommand\doubletildesigma{\skew2\tilde{\skew2\tilde{\Sigma}}}
%       italic Sigma with double tilde

% array strut to make delimiters come out right size both ends
\newsavebox{\astrutbox}
\sbox{\astrutbox}{\rule[-5pt]{0pt}{20pt}}
\newcommand{\astrut}{\usebox{\astrutbox}}

\newcommand\GaPQ{\ensuremath{G_a(P,Q)}}
\newcommand\GsPQ{\ensuremath{G_s(P,Q)}}
\newcommand\p{\ensuremath{\partial}}
\newcommand\tti{\ensuremath{\rightarrow\infty}}
\newcommand\kgd{\ensuremath{k\gamma d}}
\newcommand\shalf{\ensuremath{{\scriptstyle\frac{1}{2}}}}
\newcommand\sh{\ensuremath{^{\shalf}}}
\newcommand\smh{\ensuremath{^{-\shalf}}}
\newcommand\squart{\ensuremath{{\textstyle\frac{1}{4}}}}
\newcommand\thalf{\ensuremath{{\textstyle\frac{1}{2}}}}
\newcommand\Gat{\ensuremath{\widetilde{G_a}}}
\newcommand\ttz{\ensuremath{\rightarrow 0}}
\newcommand\ndq{\ensuremath{\frac{\mbox{$\partial$}}{\mbox{$\partial$} n_q}}}
\newcommand\sumjm{\ensuremath{\sum_{j=1}^{M}}}
\newcommand\pvi{\ensuremath{\int_0^{\infty}%
  \mskip \ifCUPmtlplainloaded -30mu\else -33mu\fi -\quad}}

\newcommand\etal{\mbox{\textit{et al.}}}
\newcommand\etc{etc.\ }
\newcommand\eg{e.g.\ }


% My Personal new commands
%%%%%%%%%%5
\def\bea{\begin{eqnarray}}
\def\eea{\end{eqnarray}}
\def\nn{\nonumber}
\def\bu{\mathbf{u}}
\def\bt{\mathbf{t}}
\def\hsigma{{\hat \sigma}}
\def\hx{{\hat x}}
\def\hy{{\hat y}}
\def\hz{{\hat z}}
\def\hs{{\hat \sigma}}
\def\hphi{{\hat \varphi}}
\def\bn{\mathbf{n}}
\def\bV{\mathbf{V}}
\def\bx{\mathbf{x}}
\def\bE{\mathbf{E}}
\def\bA{\mathbf{A}}
\def\bI{\mathbf{I}}
\def\d{\mbox{d}}
\def\dd{\mbox{d}}
\def\D{\mbox{D}}
\def\DD{\mbox{D}}
\def\Ma{Ma}
\def\Ca{Ca}
\def\Bo{Bo}
\def\stress{\mbox{\boldmath{$\sigma$}}}
\def\vort{\mbox{\boldmath{$\omega$}}}
\def\Vort{\mbox{\boldmath{$\Omega$}}}
\def\ri{\mathrm{i}}
\def\hk{\hat k}
\def\yh{\hat y}
\def\ha{\hat a}
\def\hh{\hat h}
\def\olu{\overline u}
\def\olv{\overline v}
\def\gam{\gamma}
\def\lam{\lambda}
\def\kap{\kappa}
\def\pl{\partial}
\def\hr{\hat \sig}
\def\hb{\hat b}
\def\hh{\hat h}
\def\Ma{M\!a}
\def\kh{\hat k}
\def\calU{\mathcal{U}}
\def\etal{\mbox{\it et al.}}
\def\tsigma{\tilde \sigma}
\def\ts{{\tilde \sigma}}
\def\tx{{\tilde x}}
\def\tp{{\tilde \varphi}}
\def\p{\varphi}
\def\eps{\epsilon}
\def\veps{\varepsilon}
\def\cI{\mathcal{I}}
\def\iint{\int \hspace{-0.125in} \int}
\def\cV{\mathcal{V}}
\def\cF{\mathcal{F}}
\def\cD{\mathcal{D}}
\def\cS{\mathcal{S}}
\def\gra{$V_y/U_y^\infty$}
\def\grg{$2\Omega_z/\omega_z^\infty$}
\def\fr{\frac}
\def\pr{$Pr \ $}
\def\ri{$Ri_0 \ $}
\def\rib{$Ri_B \ $}
\def\re{$Re \ $}
\def\ras{$Ra^{ust} \ $}
\def\ti{$\times$}
\def\tma{$t_{m} \ $}
\def\tmd{$t_{md} \ $}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newtheorem{lemma}{Lemma}
\newtheorem{corollary}{Corollary}

\title[Supplementary Materials]{The ``zoo'' of secondary instabilities precursory to stratified shear flow transition. Part I: shear aligned convection, pairing, and braid instabilities
\\
\huge{Supplementary Materials}}

\author[A. Mashayek, W. R. Peltier]%
{A.\ns  M\ls A\ls S\ls H\ls A\ls Y\ls E\ls K\ls $^1$%
\ls \and W.\ns R.\ns P\ls E\ls L\ls T\ls I\ls E\ls R\ls $^2$}

\affiliation{$^{1,2}$Department of Physics, University of Toronto,
Ontario,  M5S 1A7, Canada\\[\affilskip]}
%$^2$Department of Physics, University of
%Toronto, Ontario,  M5S 1A7, Canada}

\pubyear{2011}
\volume{}
\pagerange{1--1}
\date{?? and in revised form ??}
%\setcounter{page}{1}

\begin{document}
\linenumbers %Used for putting line numbers
\maketitle
\pagebreak
Throughout this supplementary document, all references to the main document will be prefixed by `$Main-$'.

\section{Application of model of Corcos and Sherman (1976)}
Figure \ref{CS76}$(a)$ shows the result of the {\it CS76} calculation for a sample test case characterized by \ri$=0.08$ for the time at which the cores
have grown to their final scales. The streamlines correspond to
the final state of the Stuart vortices and the thick black line shows the calculated braid. Panel $(b)$ of the figure shows
the evolution of the half thickness of the core (or equivalently the elevation of the tip of the braid). 
Panel $(c)$ in figure \ref{CS76} shows the evolution of the stagnation point strain rate for our test case. 
As we will be primarily interested in studying the secondary instabilities occurring
after the primary KH wave has achieved its maximum amplitude, the final value of $\gamma$ will be used in our calculations.
To investigate how the calculations compare to our numerical results, we plot the calculated
braid profiles versus the simulated profiles (both at the time of maximum $H$) in panel $(d)$ of the figure. The solid lines are from 
{\it CS76}-type calculations while the dashed lines are for the equivalent braids from the DNS simulations. The arrows in panel $(d)$ indicate the direction of increase in the
Richardson number from $0.04$ to $0.12$, $0.16$, and $0.2$ respectively. The results presented in the figure demonstrate excellent agreement between the calculations
and the simulations. As originally pointed out by {\it CS76}, the differences between the solid and dashed curves are partly due to 
simplifications made to obtain equations ($Main-3.5, \; 3.6$) and also to advecting the vorticity field using the simple
structure of the Stuart vortices. The increasing difference between the curves for larger values of the Richardson number
is derivative of the fact that we have used a value of $\Delta \rho =2$ in equation ($Main-3.6$). This value is more accurate for small \ri cases than for large. 
As \ri increases, the time of roll-up 
of the KH wave to its maximum amplitude increases and hence the density difference across the braid decreases due to the direct influence of diffusion. 
%%%%%%%%%%%%%
\begin{figure}%[h!]
\begin{center}

\vspace{0.125in}
\hspace{-2.0in} ({\it a})
\hspace{2.0in} ({\it b})

\epsfig{file=Figs/CS76/sample_streamline_Ri_08.eps,width=2.50in}
\epsfig{file=Figs/CS76/sample_case_H_Ri_08.eps,width=2.50in}
                                                                                   
\vspace{0.125in}
\hspace{-2.0in} ({\it c})
\hspace{2.0in} ({\it d})

\epsfig{file=Figs/CS76/sample_case_gamma_Ri_08.eps,width=2.50in}
\epsfig{file=Figs/CS76/braid_evolution.eps,width=2.50in}

%\begin{description}
%\item[Figure 2]
\caption{$(a)$ Streamlines and braid location obtained from {\it CS76} calculations for \ri=0.08; $(b)$ Time evolution 
of the braid elevation $H$ (or half the core thickness) for the case of panel $(a)$;
$(c)$ Time evolution of the stagnation point strain rate for the case of panel $(a)$;
$(d)$ Comparison of braids calculated using the {\it CS76} method (solid lines) and braids from simulations
at $Re=1000$ and $Pr=1$ (dashed lines). The arrows indicate increase in \ri from $0.04$ to $0.12$, $0.16$ and $0.2$.
}
%\end{description}
\label{CS76}
\end{center}
\end{figure}
%================

\section{Heuristic model for vortex core}
\subsubsection{Estimating \tma}
Figure \ref{tma}$(a)$ plots the values
of \tma versus the Richardson number for various Reynolds and Prandtl numbers (for the cases listed in $Main-table \; 1$). The overlap of
essentially all of the curves confirms the \ri-dependence of \tma. The level of stratification 
of the initial shear layer sets the maximum amplitude that the KH wave can obtain and the time it takes to reach this amplitude. 
For comparison, \tma calculated based on the {\it CS76} method 
is also plotted in the figure (shown by the $-o-$ line). As discussed earlier, evolution timescale is not accurately predicted using
their method.
%%%%%%%%%%%%%
\begin{figure}
\begin{center}

\hspace{-0.0in} ({\it a})
\hspace{2.05in} ({\it b})

\epsfig{file=Figs/t_ma/t_max.eps,width=2.55in}
\epsfig{file=Figs/del_rho_SAR/del_rho_SAR.eps,width=2.55in}


\caption{$(a)$Simulation-based results obtained for \tma
for $Pr=1, \ Re=1000$ (solid line with '+'s), $Pr=1, \ Re=2000$ (solid line with '*'s), 
$Pr=2, \ Re=1000$ (solid line with upward triangles), $Pr=2, \ Re=2000$ (solid line with downward triangles), 
$Pr=4, \ Re=1000$ (solid line with diamonds) and $Pr=4, \ Re=2000$ (solid line with crosses). The solid line with
hollow circles shows the \tma results obtained by using the {\it CS76} model;
$(b)$Simulation-based results for $\Delta \rho^{ust}_m$ versus \tmd
for $Pr=1, \ Re=1000$ (solid line with 'o's), $Pr=1, \ Re=2000$ (solid line with '+'s), 
$Pr=2, \ Re=1000$ (solid line with '*'s), $Pr=2, \ Re=2000$ (solid line with diamonds), 
$Pr=4, \ Re=1000$ (solid line with upward triangles) and $Pr=4, \ Re=2000$ (solid line with downward triangles). The solid straight line is a linear fit to all the curves in the figure. 
}
%\end{description}
\label{tma}
\end{center}
\end{figure}
%%%%%%%%%%%%%

\subsubsection{Estimating $\Delta \rho^{ust}_m$}
The results listed in $Main-table \;  1$ for $\Delta \rho^{ust}_m$ are plotted in figure \ref{tma}$(b)$. The markers on each of the lines in the figure correspond to $Ri(0)=0.04, \ 0.08, \ 0.12, \ 0.16,$ and $0.2$. The figure shows that $\Delta \rho^{ust}_m$ is very close to $2$ for small \ri and drops to a value close to $0.5$ at $Ri(0)=0.2$. It also shows that for higher Reynolds and Prandtl numbers, $\Delta \rho^{ust}_m$ decreases less significantly with \ri. So, for very high values of \re and \pr (outside the range covered by the cases explicitly analyzed in our study) it might be safe to assume a value of $2$ for $\Delta \rho^{ust}_m$. In our calculations of the \ras however, we will use the linear fit shown in the figure, the equation for which is:
\bea
\Delta \rho^{ust}_m=3.22-0.024 t_{m}.
\eea

It should be noted that the abscissa in figure \ref{tma}$(b)$ is labled $t_{md}$ which refers to the time at which $\Delta \rho^{ust}$ peaks. Although table $SM-1$ shows that \tma and \tmd are close, it is important to note that they are not exactly the same. Figure \ref{t_ma_t_md}$(a)$ shows the ratio of these two times and indicates that for $0.04<Ri(0)<0.16$, it is reasonable to use \tma for \tmd. This will be our choice as a simplifying step in the analysis to follow. Comparison between the \tma and \tmd columns in table $Main-table \;  1$ shows that the latter is more sensitive to \pr and \re. This is because the time needed by the KH wave to grow to its maximum amplitude is primarily a function of the Richardson number while evolution of the unstable region is affected by the inter-layer diffusion inside the core which itself depends on both \pr and \re.
%%%%%%%%%%%%%
\begin{figure}
\begin{center}

\hspace{-0.0in} ({\it a})
\hspace{2.05in} ({\it b})

\epsfig{file=Figs/del_rho_SAR/t_ma_over_t_md.eps,width=2.55in}
\epsfig{file=Figs/del_SAR/del_SAR.eps,width=2.55in}

\caption{$(a)$ $t_{m}/t_{md}$ plotted for the cases listed in a table in $Main-table \; 1$. At each value of \ri, there are six stars corresponding to different \re and \pr; $(b)$Data points from the same table (circles) plotted along with correlation ($Main-(3.3)$). The Richardson number increases from $0.04$ for the data points at the left to $0.2$ for points at the right.
}
%\end{description}
\label{t_ma_t_md}
\end{center}
\end{figure}
%%%%%%%%%%%%%


\section{Analysis of the braid}

Table \ref{t2} shows the results of the theoretical prediction of the braid stagnation point Richardson number versus those obtained from 2D numerical simulations of the table in $Main-table \; 1$. 
To obtain theoretical values for \rib, we have employed ($Main-3.9$) along with $\gamma$ and $\psi$
obtained using the {\it CS76} model (using equations ($Main-3.7, \; 3.8$). 
It is important to note that in the analyses of {\it CS76} leading to relation $SM-(2.1)$, it was assumed that $\gamma$ could be assumed to be constant along the braid. However, the main equation governing evolution of $\gamma$ in {\it CS76} indicates a non-constant value with a peak at the stagnation point. In our analysis, we average $\gamma$ over the braid and use the result in ($Main-3.9$). 
For comparison between the predictions of ($Main-3.9$) and the simulations, we choose the braid Richardson number at the stagnation point ($Ri_B^{stag}$) in our numerical results and at a time between \tma and the onset of the pairing instability. As we will see, $Ri_B^{stag}$ does not change significantly during this period of time. 
Table \ref{t2} shows that the values of $Ri_B^{stag}$ calculated using the theoretical prediction are in good agreement with our 2D numerical simulations. The agreement however diminishes with an increase in the Prandtl number. 
This is probably because, for higher Prandtl numbers, the strength of the vorticity bands in the cores is increased compared to that for the lower Prandtl number cases. Hence, the effect of the flow field of the cores on the braid at the stagnation point is more pronounced for cases with $Pr=2$. As the theoretical prediction totally ignores this effect, the agreement between the theory and the simulations deteriorates at higher Prandtl numbers. Nevertheless, the results are still in reasonable agreement for $Pr=2$.

\begin{table} 
\begin{center}
\begin{tabular}{l|cccccccccccccccccccccccccccccccccccc}
\hline
Case Name       &&& $Ri_B$ theory &&& $Ri_B^{stag}$ simulation\\
$c1-1000-0.08$  &&& 0.20  	    &&& 0.20\\
$c1-1000-0.16$  &&& 0.12  	    &&& 0.13\\
$c1-2000-0.08$  &&& 0.14  	    &&& 0.15\\
$c1-2000-0.16$  &&& 0.085 	   &&& 0.09\\
$c1-4000-0.08$  &&& 0.10  	   &&& 0.13\\
$c1-4000-0.16$  &&& 0.06  	   &&& 0.04\\
$c2-1000-0.08$  &&& 0.36  	   &&& 0.25\\
$c2-1000-0.16$  &&& 0.21  	   &&& 0.15\\
$c2-2000-0.08$  &&& 0.25  	   &&& 0.19\\
$c2-2000-0.16$  &&& 0.14  	   &&& 0.12\\

\hline
\end{tabular}
\caption{\rib calculated using equation ($Main-3.7$) in the middle column and from simulations in the right column for sample test cases from $Main-table \; 1$}
\label{t2}
\end{center}
\end{table}





\subsection{Diagnostic tools for the braid}
To facilitate our analyses of the braid, we needed to develop certain diagnostic tools. To explain development of the tools, we consider the case c1-2000-0.12 as an example. Figure \ref{braid_intro} shows the vorticity field at a sequence of times in
the evolution of this flow.
%Panel $(d)$ shows that once the intense vorticity bands of the core extend to become sufficiently close to the braid, an instability is excited (frame $(e)$). This instability, which will be discussed in more details later, leads to deformation of the braid and to roll-up of two small secondary vortices
%on it (frame $f$). At the time that the cores' vorticity bands reach the braid, the vorticity in the braid is drained as shown in panels $(f-j)$. 
%The existence of the braid plays an important role in maintaining  the large scale structure
%of the KH wave in a stratified layer.\footnote{As explained in {\it CS76}, the height of vortex cores in a KH billow saturate once the vorticity in the braid is transported to the cores and there is a balance between baroclinic production of vorticity in the braid and loss of vorticity in the cores. This is not the case in the unstratified case.} In a stratified layer, once this vorticity band is drained and the strain rate along the braid has decreased significantly,
%the two cores begin to orbit each other as two isolated vortices with the same sense of rotation would do when they are sufficiently proximate.
%This marks the beginning of the pairing instability (panels $(j-l)$). Once the pairing process begins, a 
%new braid is formed once more (panel $(l)$) between the large amplitude cores that make up the train of billows. 
%The tilt angle of the braid after the two cores have merged (panel $(l)$)
%is close to that of the initial fully-developed KH wave (panel $(b)$) because both the scale and the distance between the cores doubles during the pairing process. The tilt angle
%of the braid is governed by \ri both before the pairing instability begins and after the two cores have completely merged.
%\begin{landscape}
%%%%%%%%%%%%%
\begin{figure}
\begin{center}

\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_50.eps,width=2.56in}
\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_60.eps,width=2.56in}

\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_75.eps,width=2.56in}
\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_80.eps,width=2.56in}

\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_85.eps,width=2.56in}
\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_90.eps,width=2.56in}

\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_95.eps,width=2.56in}
\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_100.eps,width=2.56in}

\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_105.eps,width=2.56in}
\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_110.eps,width=2.56in}

\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_140.eps,width=2.56in}
\epsfig{file=Figs/braid/rho_omeg_plots/omeg_t_165.eps,width=2.56in}
\caption{Vorticity contours for case c1-2000-0.12 and at times $(a)$ t=50, $(b)$ t=60, $(c)$ t=75, $(d)$ t=80, $(e)$ t=85, $(f)$ t=90, $(g)$ t=95, $(h)$ t=100, $(i)$ t=105, $(j)$ t=110, $(k)$ t=140, $(l)$ t=165.  
}
%\end{description}
\label{braid_intro}
\end{center}
\end{figure}
%================
%\end{landscape}

To extract braid information from the simulations, we extracted data from one hundred traverses between the centers of the cores for each time frame of each simulation. The point of maximum 
vorticity on each traverse was then identified to define the location of the braid center. Alternatively, one could use the density profiles to locate the center of the braid. Connecting the $100$ points obtained in this fashion provided a braid
profile for each time frame. The solid lines in Figure \ref{braid_diag}$(a)$ present the results for three different times for the same test case as illustrated in Figure \ref{braid_intro}.
The solid lines in the figure are overlain by a third-order polynominal fit shown as the dash-dotted curve. For times shown in panel $(a)$, the braid profiles and the fits coincide closely. Panel $(b)$ illustrates the braid at a sequence of later times when the braid deformation is pronounced and so the fits do not coincide with the earlier braid profiles. For these cases, the fits (dash-dotted lines) can be treated as an unperturbed braid, or in other words, as the braid had it had no secondary braid instabilities deforming it.  Therefore, for the cases of panel $(b)$, the deviation of the braids from the polynominal fits can be considered a measure of the braid deformation. We have calculated the area between the solid and dash-dotted curves for each time step and we employ the difference to define the ``braid deformation''. For the test case under consideration, this braid deformation is plotted in panel $(c)$ of Figure \ref{braid_intro}. 
At early times when the braid has yet to become deformed, this parameter has a value near zero. However, at a time shortly before $t=80$, the braid begins to grow (as shown in Figure \ref{braid_intro} $(d-e)$). At a time slightly prior to $t=90$, small billows form on the braid and these appear in the form of small peaks on the braid deformation plot. The larger deformations at later times in panel $(c)$ of figure \ref{braid_diag} are associated with the pairing process. Comparing panels $(a,b)$ shows that despite the braid deformation, its location does not vary much for times beyond $t=50$.


Panel $(d)$ of Figure \ref{braid_intro} shows the braid thickness calculated from the vorticity field (solid line) and the density field (dashed line). We define the braid vorticity thickness as twice the normal-to-braid distance of the center of the braid from the location where the vorticity has a value of about $10 \%$ of the difference between the maximum value (at the center of the braid) and the minimum value of approximately zero (sufficiently far from the braid). The density thickness is defined in a similar fashion. This definition of the thickness based upon the $10 \%$ criteria is somewhat arbitrary and a more accurate definition is that proposed in {\it S03}. However, we are primarily interested in braid location and changes in thickness not the value of braid thickness itself. As shown in the figure, there is a period of time between $t=50$ and $t=85$ when the braid is in a stable condition and its thickness is unchanging.
As we will see, the value of \rib is also relatively stable during this period and, if it is sufficiently low, we may anticipate the emergence of secondary shear instability of the braid. The ratio of the braid's vorticity thickness to its density thickness is also close to unity during the semi-equilibrium phase which is in agreement with the results of {\it S03} for $Pr=1$. 


Panel $(e)$ plots the parameter $\delta_{pair}$ which is the horizontal distance between the tips of the strong vorticity layers in two neighbour cores (shown in panels $(b-d)$ of Figure \ref{braid_intro}). As mentioned previously, the time at which these bands meet and thereafter drain the vorticity from the braid is of critical importance as it marks the onset of the growth of various secondary instabilities. As panel $(e)$ of Figure \ref{braid_diag} shows, at a time shortly after $t=80$, $\delta_{pair}$ drops to zero and it is at this time that we observe the braid deformation to begin to grow in panel $(c)$. 
This time also marks the onset of the merging instability. To demonstrate this, we plot $\Delta z_{core}$ in panel $(f)$ which is the difference between the vertical positions of tops of the left and right cores. Prior to the onset of merging, the two positions are identical and $\Delta z_{core}$ is zero but this grows slowly after $\delta_{pair}$ vanishes implying that pairing is underway.
Although the coincidence of the time of $\delta_{pair}=0$ and the onset of $\Delta z_{core}>0$ occurs in many of our simulations, in general, the pairing instability may get underway even before $\delta_{pair}$ vanishes depending on the flow initialization. Moreover, $\delta_{pair}=0$ is not a necessary condition for the merging instability. At high Richardson numbers ($0.2$ and higher) merging may be completely prohibited by the influence of the density stratification. 


Panel $(g)$ of Figure \ref{braid_intro} presents time variation of the vertical extent of the core, $D_{core}$, demonstrating that it oscillates about the mean for a time but is thereafter affected by the onset of secondary instabilities that originate on the braid and by the pairing process. Panel $(h)$ plots the variation of \rib at the stagnation point. During the period when the braid is approximately steady ($t \sim 50-85$), the $Ri_{B}^{stag}$ is essentially constant. Near $t=85$ when $\delta_{pair}$ vanishes and the braid is drained of its vorticity, its Richardson number drops very rapidly and secondary vortices appear on the braid.
The large variations for $t>100$ in panel $(h)$ are due to the merging process. 
A more detailed examination of \rib may be constructed by plotting it as a function of position along the braid.
This is done in panels $(i)$ and $(j)$ for various times. For earlier times and before $\delta_{pair}$ tends to zero, the minimum value of \rib obtains at the stagnation point. However, as $\delta_{pair} \rightarrow 0$, at the locations where the tips of the cores outermost negative vorticity bands meet the braid, \rib is further reduced. This is clearly seen for times $t \ge 75$. The two points of contact have a near zero Richardson number at $t=85$ which is near the time that secondary vortexes form at those points on the braid. The horizontal axes in panels $(i)$ and $(j)$ denote the global x-direction and not the along-braid coordinate. 
The diagnostic tools described above (specially those employed to construct panels $(c)$ through $(f)$) are useful for determining conditions under which the emergence of small and large scale deformations on the braid and cores appear. We will further employ these tools in our analyses in the next sections and we will also construct similar plots to those shown in figure \ref{braid_diag} for the cases to be considered.

%\begin{landscape}
%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\hspace{-0.3in} ({\it a})
\hspace{1.6in} ({\it b})

\epsfig{file=Figs/braid/braid_location_ri/braid_location1.eps,width=2in}
\epsfig{file=Figs/braid/braid_location_ri/braid_location2.eps,width=2in}

\hspace{-0.3in} ({\it c})
\hspace{1.6in} ({\it d})

\epsfig{file=Figs/braid/braid_deform.eps,width=2in}
\epsfig{file=Figs/braid/braid_thickness_1_third.eps,width=2in}

\hspace{-0.3in} ({\it e})
\hspace{1.6in} ({\it f})

\epsfig{file=Figs/braid/Del_omeg.eps,width=2in}
\epsfig{file=Figs/braid/del_y_cores.eps,width=2in}

\hspace{-0.3in} ({\it g})
\hspace{1.6in} ({\it h})

\epsfig{file=Figs/braid/D_left_core.eps,width=2in}
\epsfig{file=Figs/braid/Ri_stag.eps,width=2in}

\vspace{0.25in}

\hspace{-0.3in} ({\it i})
\hspace{1.6in} ({\it j})
                                                                                  
\epsfig{file=Figs/braid/braid_location_ri/braid_Ri1.eps,width=2in}
\epsfig{file=Figs/braid/braid_location_ri/braid_Ri2.eps,width=2in}

%\begin{description}
%\item[Figure 2]
\caption{Flow diagnostic tools for case c1-2000-0.12. $(a)$ Braid location at times $t=30, \ 50, \ 75$; $(b)$ braid location at times $t=80, \ 85$;
$(c)$ braid deformation; $(d)$ braid thickness calculated from the vorticity and the density fields; $(e)$ $\delta_{pair}$; $(f)$ $\Delta z_{cores}$; 
$(g)$ $D_{core}$; $(h)$ $Ri_B$ at the stagnation point; $(i)$ $Ri_B$ along the braid for $t=30, 50, \ 75$; $(j)$ $Ri_B$ along the braid for $t=80, 85$.}

%\end{description}
\label{braid_diag}
\end{center}
\end{figure}
%================
%\end{landscape}

\input{Supplement_Stability_analysis}	 		 
%\pagebreak
\bibliographystyle{jfm}
% Note the spaces between the initials
\bibliography{Mashayek_Peltier_2012_paper1}
\end{document}
