Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Linear stability analysis of a oceanic frontal system over a Submarine ridge De︠ry, Francis 1997

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1997-0227.pdf [ 5MB ]
JSON: 831-1.0079905.json
JSON-LD: 831-1.0079905-ld.json
RDF/XML (Pretty): 831-1.0079905-rdf.xml
RDF/JSON: 831-1.0079905-rdf.json
Turtle: 831-1.0079905-turtle.txt
N-Triples: 831-1.0079905-rdf-ntriples.txt
Original Record: 831-1.0079905-source.json
Full Text

Full Text

Linear Stability Analysis of a Oceanic Frontal System over a Submarine Ridge By Francis Dery B.Sc. of Physics, Universite de Montreal M.Sc. of Atmospheric and Oceanic Sciences, McGill University A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER 6f SCIENCE -in THE FACULTY OF GRADUATE STUDIES Department of MATHEMATICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April 23, 1997 ©Francis Dery, 1997 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Mathematics The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1Z1 Date: ipaLLh. ABSTRACT The stability of frontal ly-trapped baroclinic w a v e s along the crest of an oceanic ridge on an f - p l a n e is examined to determine under which condit ions two neighbor ing, inviscid water layers, each of uniform density and m o m e n t u m , can coexist without mixing by eddies associated with instabilities. A l inear stability analysis for a w e a k disturbance is per formed over a s teady laminar state. The study fol lows that of Gawark iewicz [1991] for a frontal sys tem over a shelf-break, but deviates by building an implicit dispersion relation to solve, by finite dif ference, a four th -o rder ordinary differential equat ion. The implicit d ispersion relation for a submar ine ridge, with straight s lopes that ex tend to infinity, depends on four parameters that character ize the physical sys tem. They are the nondimensional w a v e number , the nondimensional ized ridge s lopes, and a quantity, ana logous to the Richardson number , that character izes the shear ing force against the buoyancy at the pycnocl ine. For a model with finite straight s lopes, the posit ions at which the s lopes meet the flat bot tom have to be taken into account in the boundary condit ions. The mode l is appl ied to var ious cases. First, both s lopes are set to zero for a flat topo-graphy case and for compar ison with Orlanski [1968] . In the second case, a s lope is introduced beneath the front to s imulate a shelf-break analogous to Flagg and Beardsley 's [1978] model with a seaward false bot tom condit ion. Stable branches of backward -p ropa -gat ing w a v e s occur with phase velocity faster than the layers' speed . Unstable features occur for w a v e s with phase velocity less or equal to the layers' speed . In further cases, the false bot tom condit ion is replaced by shifting the foot of the ridge away f rom the shelf break, into a near or far location. The main f indings are that a r idgefoot d isp lacement does not change most unstable features, but it increases the number of stable branches. Finally, the shel fbreak model is modif ied for a symmetr ic r idge mode l by introducing another finite straight slope of the s a m e value and length. In these cases, the main unstable features are not different f rom the shel fbreak model . It creates fo rward-propagat ing w a v e s with phase velocity faster than the layers' speed . iii TABLE OF CONTENTS Abstract » Table of contents iii List of figures iv List of tables vii INTRODUCTION 1 CHAPTER 1: Linear stability analysis 5 1.1 Formulation of the model 5 1.2 Triangular Ridge 12 1.3 The dispersion relation 16 1.4 Initial parameters 22 CHAPTER 2: Numerical solutions of the dispersion equations 24 CHAPTER 3: Results 27 3.1 Flat cases 27 3.1.1 Low Richardson number 27 3.1.2 High Richardson number 30 3.2 Flagg and Beardsley shelf-break model 33 3.3 Shelf-break model with foot at y=2 42 3.4 Symmetric slopes ridge with feet at y=-1 and y=2 49 3.5 Shelfbreaks and symmetric slopes ridges with long slope 55 3.5.1 Shelfbreak with foot at y=10 and R/=57 55 3.5.2 Symmetric slopes ridge with R/=57, feet at y=-9 and y=10 60 CHAPTER 4: The Lomonosov Ridge 63 4.1 The observed front during ODEN 91 expedition 63 4.2 The modeled ridge with a blunt crest 67 4.3 The effect of a sharper crest on the Lomonosov Ridge 72 CHAPTER 5: Discussion and conclusion 75 Bibliography 79 APPENDIX A: The Tricomi solution of the confluent hypergeometric equation 81 APPENDIX B: Subroutines to compute the dispersion relation and search of roots 89 iv LIST OF FIGURES Figure 1.1 .s 6 Two layers bordered by a submarine ridge, with density px,p2, (the denser layer in region III), alongridge velocities UX,U2, and thicknesses r\{y), h^iy), respectively. The top of the ridge is at a depth H. The frontal width is L and depends on velocity shear and density difference. The thick oblique line starting from the crest represents the pycnocline in the laminar situation. Figure 1.2 12 Same as fig. 1.1, except that layer thicknesses are defined as }\ {y) = H(1 - yly) in region I, as K 00 = H(l + y^y) in region III, and, in region II, as i\ 00 = H(l -y), r\ (y) = H{\ + y2 )y for the upper and lower layers respectively, STJ is the displacement of the perturbed pycnocline. Figure 3.1 28 Flat bottom and Ri=3. The upper graph represents a nondimensional Doppler-shifted frequency against wave number. Al l modes are contained inside a range (~Ro,Ro) shown by dashed lines. The lower graph is nondimensional growth rate. Figure 3.2 29 Flat bottom and Ri = 3. The Orlanski's eigenvalue representation with Ro as a state variable. Due to symmetry, only the positive side of the real part is shown. Figure 3.3 31 Flat bottom configuration with /?/=57. On the left panel: the observed stable points (up) and their x-transformation (down). On the right panel: the observed unstable points with frequency (up) and growth rate (down). Due to symmetry, only positive real parts are shown. Figure 3.4 32 The x-representation at Ri = 57 for the flat bottom configuration. Figure 3.5 34 FB configuration with y2= 0.25. Left panel: the frequency of stable modes outside the Rossby range (up) and the x-transform of stable branches inside the Rossby range (down). The dashed line is the lower Rossby barrier. Right upper panel: the real part of the unstable points. Dashed lines are the upper and lower Rossby barriers. Right lower panel: the imaginary part of the unstable points. Figure 3.6 35 Same as fig. 3.5, but for a FB model with y2=0A0. Figure 3.7 36 Same as fig. 3.5, but for a FB model with /2=0.63. V Figure 3.8 37 Same as fig. 3.5, but for a FB model with y2=\.0. Figure 3.9 38 Same as fig. 3.5, but for a FB model with y2=l.6. Figure 3.10 39 Same as fig. 3.5, but for a FB model with y2=2.5. Figure 3.11 40 Critical wave number for parabolic stable branches Kc vs y2. Figure 3.12 42 Same as fig. 3.5, but for a shelfbreak with /2=0.25 and foot at y=2. Figure 3.13 43 Same as fig. 3.5, but for a shelfbreak with y2=QA and foot at y=2. Figure 3.14 44 Same as fig. 3.5, but for a shelfbreak with y2=0.63 and foot at y=2. Figure 3.15 45 Same as fig. 3.5, but for a shelfbreak with ^2=1.0 foot at y=2. Figure 3.16 46 Same as fig. 3.5, but for a shelfbreak with y2=l.6 and foot at y=2. Figure 3.17 47 Same as fig. 3.5, but for a shelfbreak with y2=2.5 and foot at y=2. Figure 3.18 48 Upper: Interpolated group velocities for a very long wave. Lower: for the second type of stable mode. The absciss is the shelfbreak slope y2. Note that the scales for the two curves are different in the upper panel. Figure 3.19 49 Sketch of the topography for a ridge of symmetric slopes and feet respectively at y=-l and y=2 Figure 3.20 50 Same as fig. 3.5, but for a symmetric slopes ridge with y=0.25 and feet at y=-l and y=2. Figure 3.21 51 Same as fig. 3.5, but for a symmetric slopes ridge with v=0.40 and feet at y=-l and y=2. Figure 3.22 52 Same as fig. 3.5, but for a symmetric slopes ridge with y=0.63 and feet at y=-l and y=2. vi Figure 3.23 52 Same as fig. 3.5, but for a symmetric slopes ridge with 7=1.0 and feet at y=-l and y=2. Figure 3.24 53 Same as fig. 3.5, but for a symmetric slopes ridge with 7=1.6 and feet at y=-l and y=2. Figure 3.25 54 Same as fig. 3.5, but for a symmetric slopes ridge with 7=2.5 and feet at y=-l and y=2. Figure 3.26 56 Same as fig. 3.5, but for a shelfbreak with /2=0.25 and foot at y=10. Figure 3.27 57 Same as fig. 3.5, but for a shelfbreak with /2=0.4 and foot at y=10. Figure 3.28 58 Same as fig. 3.5, but for a shelfbreak with /2=1.0 and foot at y=10. Figure 3.29 59 Same as fig. 3.5, but for a shelfbreak with /2=2.5 and foot at y=10. Figure 3.30 61 Symmetric ridge with Ri=51, 7=1.0 and feet at y=10 and y=-9. Upper panel: the stable modes outside the Rossby region. Lower panel: the x-transformation of the frequency for stable points inside the Rossby region. Figure 3.31 62 Symmetric ridge 7=1.0 and feet at y=10 and y=-9. Right upper panel: the real part of the unstable points. Dashed lines are the upper and lower Rossby barriers. Right lower panel: the imaginary part of the unstable points. Figure 4.1 64 Contour plots of potential temperature observed during the ODEN 91. The dashed line from the top of the Lomonosov Ridge shows the location of the front beneath the pycnocline Figure 4.2 65 Contour plots of potential temperature observed during the ODEN 91 in the top 1000 m. The dashed lines show vortex sheets in the neighborhood of the Lomonosov Ridge. The vortex sheet on the right is the frontal interface beneath the pycnocline. The arrows show a sinking process of water of Atlantic origin into the Canadian water assembly. VII Figure 4.3 65 Profile of potential temperature (T), salinity (S) and density-1000 (a) between depths of 220 decibars and 320 decibars from station 23 of ODEN 91. Arrows indicate levels where vertical gradient of density is positive, which triggers vertical mixing. The water column between 236 and 262 decibars is going to undergo vertical mixing. Figure 4.4 68 The model representation of the Lomonosov Ridge frontal system. Figure 4.5 69 The real and imaginary parts of the eigenvalue for the unstable mode in the Lomonosov Ridge case with flattened crest. Figure 4.6 70 The stable modes outside the Rossby range for the Lomonosov Ridge with a flattened crest. Figure 4.7 71 The stable eigenvalues with the x-transformation inside the Rossby range for the Lomonosov Ridge with flattened crest. The regions 0.75<|x|<l were avoided because too compact for interpretation. Figure 4.8 72 Stable modes outside the Rossby range for a narrow Lomonosov Ridge. Figure 4.9 73 Stable modes with the x-transformation inside the Rossby range for a sharp Lomonosov Ridge. Figure 4.10 74 The unstable modes found for the sharp crest model of the Lomonosov Ridge: a) the Doppler-shifted frequency divided by the Rossby number, b) the growth rate. LIST OF TABLES Table 1 41 Group velocities of the fastest mode of first type and of the second type mode interpolated for very long waves w.r.t. the slope. The intercepts' separation is the extrapolated value of the frequency for the second mode for long waves. Table 2 Same as table 1, but the shelf-break foot is displaced to y=2. 47 1 INTRODUCTION Fronts are frequent synoptic features of marine environments that may hinder horizontal transfers of heat, salt, momentum and other properties like pollutants. The time scale over which a front is stable plays a crucial role in oceanic circulation and mixing processes. When instability occurs, the disturbance grows and distorts the front. Resulting intrusions enhance oceanic dissipation by diffusion. The front can break with the formation of eddies and rings that add to crossfrontal exchanges. Comprehension of frontal systems and their stability is needed especially for regions that are poorly monitored. One example apparent from the sparse observations under the Arctic polar icepack (Anderson et al. [1994]) is a front dividing the water assemblies of the Canadian Basin from the Eurasian Basin, about the Lomonosov Ridge. The survey of the various pollutants leaking from the former Soviet junkyards and the worries about sunken nuclear waste into the Arctic Sea requires knowledge of whether the observed front is a temporary or a permanent feature. This work has been motivated by a supposition that the frontal system can be stabilized by the Lomonosov Ridge like shelf-break fronts in more monitored areas. Recent observations from the Larsen-93 expedition reported characteristics of the Eurasian Basin water assembly below 2000 m in the Makarov basin, which is a part of the Canadian Basin (McLaughlin et al. [1996]). This was interpreted as an eastward displacement of the front on the Eurasian side, into the Canadian Basin to at least the deep layer, over the neighboring Mendeleyev-Alpha Ridge. The Larsen-93 data support a scheme of communication between the Canadian and Eurasian basins through the limits of the Lomonosov Ridge proposed by Anderson et al. [1994]. Rudels et al. [1994] proposed that water of the Eurasian Basin entering into the Canadian Basin north of the Laptev Sea is modified by colder and fresher shelfwater sinking and intruding in the transition zone beneath the Atlantic layer. The modified Atlantic water continues until the Mendeleyev Ridge, where one branch of the assembly is redirected poleward. McLaughlin et al. [1996] concludes that Larsen-93 data imply the front is located near the Mendeleyev-Alpha ridge system and notes that the Arctic marine 2 front seems to favor a position over a ridge topography. This answers partly our interrogations about stability. However, we still do not understand fully how fast a front can break by baroclinic instability coupled to its underlying topography. Frontal systems are an important aspect of geophysical fluid dynamics, both in the atmo-sphere and in the ocean. Investigations of them have looked into separate phenomena such as baroclinic instability (Charney [1947], Eady [1949]); conversion of energy between mean and eddy flows (Fjortoft [1951]); transformation of energy by baroclinic waves in a two-layer model (Phillips [1954]); and the role of the frontal instability in the formation and development of atmospheric fronts (Bjerknes [1919], Solberg [1928], Kotschin [1932], Orlanski [1968]). Kotschin [1932] considered two incompressible, homogeneous fluids, with a shear and a slight density difference (classical Margules front type), bounded above and below by two rigid, hori-zontal planes, and found a stability criterion relating the wave number with the frontal slope. However, as the density difference decreases, the interface becomes more vertical. At the limit, the front becomes barotropic, i.e., sustained by the velocity shear only, and instability is expected, but this was not admitted by Kotschin's stability criterion. This leads Orlanski [1968] to search in Kotschin's [1932] model, for a second kind of frontal instability that depends on the Rossby and Richardson numbers. For the barotropic case, he found two stable modes of magni-tude equal to the Coriolis parameter (inertia waves), a vanishing and an unstable mode without phase velocity and with a decay/growth rate equal to the Coriolis parameter. The Kotschin's [1932] model applies well to a marine case approximated by a reduced-gravity front extending from a flat bottom to a rigid lid with layers of uniform momenta on an /-plane. The primitive equations of geophysical fluid dynamics have been generally simplified. The most popular simplified forms have been the barotropic vorticity equation, the quasigeostrophic equations and the more complex balance equations. Eliassen [1949] derived a system of hydrostatic, filtered momentum equations with the Boussinesq approximation, which neglects the substantial derivative of the ageostrophic part of the horizontal wind (or horizontal momentum in the general case). Since the total horizontal wind is retained in the advection 3 terms of the substantial derivative, the system is referred as the semigeostrophic equations by Hoskins [1975], who found them useful in the analysis of ageostrophic phenomena such as frontogenesis and nonlinear baroclinic instability. Hoskins [1976] used them to study the development of the Eady wave instability mode of a zonal flow with uniform vertical shear. Duffy [1976] applied the semigeostrophic equations to Kotschin's [1932] model and compared his results with Orlanski [1968]. He found, for the barotropic case at the limit where the Rossby number vanishes, two modes of phase velocity zero and growth/decay equal to the Coriolis parameter. He missed the stable inertia waves found by Orlanski [1968]. He saw that for Ro< 0.4, the semigeostrophic equations give conventional baroclinic instability and large-scale shear instability, although he consistently underestimated the growth rate and missed small-scale Rayleigh and Kelvin-Helmholtz instabilities. His results had been foretold by Fjortoft [1962] who saw the semigeostrophic equations filter out high-frequency phenomena. Pedlosky [1964] found, as a necessary condition for instability, for a quasi-geostrophic two-layer system, that the gradients of vorticity (or thickness) of each layer must be of opposite signs, which is necessarily the case when the interface and the bottom are of opposing slopes. Orlanski [1969] showed that, in a nearly flat bottom case, an increase in the lower layer thickness (by increasing the bottom slope) causes a decrease in the growth rate. This suggests that a finite (but small enough) sloping bottom tends to stabilize the front. Flagg and Beardsley [1978] added to Orlanski's [1968] model a realistic steep bottom slope beneath the interface, while the bottom stays flat outside the frontal region, in order to simulate a shelf-break region. They found unstable modes, for which, an increasing bottom slope decreases their growth rates, hence frontal wave modes are not uniquely generated by local baroclinic instability. They also investigated the case of topographic shear waves with a seaward, exponentially-shaped bottom. Steepness was governed by a nondimensional parameter, which leads to an eigenvalue equation that was solved numerically. In the limit where steepness disappears, their eigenvalue equation is reduced to a quartic with solutions found by Orlanski [1968]. 4 Gawarkiewicz [1991] extended the constant linear bottom slope of Flagg and Beardsley [1978] model over the seaward region and solved the boundary conditions by a confluent hypergeometric function of the second kind, following the technique of Duxbury [1963], and considered the case of topographic shear waves. He used Hoskin's [1975] quasigeostrophic momentum approximation for the perturbation (\co 2\«1) and a theorem from Barth [1989] to reduce the search of the unstable eigenvalues inside a semicircle in the complex plane, of radius equivalent to the Rossby radius. He found that barotropic shear waves were very sensitive to the choice of bottom topography and his stable waves differ from the case studied by Flagg and Beardsley [1978]. Also, he found that purely barotropic instabilities are hindered by the seaward sloping bottom topography, and that, as shown by Flagg and Beardsley [1978], increasing the bottom slope reduces the instability substantially. He found for the baroclinic case that a linearly sloping, unbounded topography in the seaward region does not substantially change the result of the analysis versus the flat seaward bottom case. The framework of this thesis will be divided as followed: In the next chapter, we will formulate the model for a general bottom topography. Then the topography will be defined as a wedge-shaped ridge. We will use a series expansion to bring the linear set of four ordinary differential equations into an implicit dispersion relation for the eigenvalue that contains the growth rate in its imaginary part. The last part of the chapter describes the free parameters of the model. The second chapter describes the schemes used to solve the dispersion relation numerically. The third chapter shows the results of the analysis applied to various cases: flat horizontal bottom model at low and high Richardson numbers; shelfbreak models with seaward flat bottom conditions or with near or remote slope foot location; symmetric slopes ridge models with near or remote locations of slope feet. The fourth chapter describes the application of the model to the Lomonosov Ridge front observed during the ODEN 91, but with a simplified water density structure. The fifth chapter will contain the final discussion. 5 CHAPTER 1 LINEAR STABILITY ANALYSIS 1.1 Formulation of the model The case that we consider is that of a velocity discontinuity along a submarine ridge between two layers with different, uniformly distributed densities and no crossfrontal velocity. The front, in the basic steady laminar state, will stretch linearly from the ridge toward the heavier side in order that pressure gradient forces counteract the Coriolis forces (see fig. 1.1). The basic properties in the upper layer, referred to as layer 1, are density, p„ alongridge velocity, and thickness, hfyr), while those in the lower layer, referred to as layer 2, are density, p2, alongridge velocity, U2, and thickness, h2(y). The crest of the submarine ridge is at a depth H under the surface and the front surfaces at a distance L from the ridge. The system is rotating at a period 4 ; r / / , where / is the Coriolis parameter, assumed constant (/"-plane approximation). We assume also there is no surface wave (rigid lid approximation). In the basic state, the pressure gradient force, exerted horizontally on a water parcel, is in geostrophic balance with the crossfrontal Coriolis force in both layers: dP. a.) — where j= 1,2 are indices of the layers. The coordinate system is defined such that the origin of the crossfrontal coordinate,^, is at the crest of the ridge, and that of the depth coordinate z is at the level where the front surfaces. The JC-direction is into the page. Figure 1.1. Two layers bordered by a submarine ridge, with density px,p2, (the denser layer in region III), alongridge velocities U1,U2> and thicknesses )\ (y), r\ (y), respectively. The top of the ridge is at a depth H. The frontal width is L and depends on velocity shear and density difference. The pycnocline in the laminar situation is represented by the dashed oblique line starting from the crest. Defining the surface pressure to be zero as a boundary condition, the integration of (1.1) yields (1.2a) Px (y,z) = pxgz - fpxUx (y-L) in layer 1, (1.2b) P2 (y, z) - p2gz - fp2U2 (y - L) in region III of layer 2. The slope of the pycnocline depends on the density contrast and on the velocity shear, which are both constant for vertically uniform layers. Thus the depth of the pycnocline relative to the level where the front crops out is H(I-y/L). The geostrophic term, which couples the layer momentum and the Earth rotation, causes the sea surface to be slightly slanted, i.e. not quite horizontal (zSUKrace = (y-L)Ujf/g). However, the order of the surface variation is negligible relative to the front and topographic feature. The pressure equation under the pycnocline (layer 2 in region II) is P2(y,z) = Pl{y,H(l-y/L)) + p2g{z - H(l-y/LJ) = pxgH{\-y/L) - fpx Ux(y-L) + p2g(z- H{\ - y/L)). The geostrophic balance under the pycnocline (layer 2, region II) yields fp2u2 = -4 t=(a -P 2 )^r+fPiU l or ay L (1.3) f(P2U2-plU1) = (p1~p2)^j-. The left-hand side of (1.3) is the product of the Coriolis parameter and the momentum shear. Using the Boussinesq approximation, the density difference contributes negligibly to the momentum shear and the geostrophic balance can be rewritten / p A < y = ( p 2 - A ) ^ , where AU=Ul-U2 is the velocity shear and p is a reference density. The representation ,of our model implies UX>U2. The frontal length can be expressed as a function of the density difference and the velocity shear / A C / ' where g' = {p2-px )g/p is the reduced-gravity parameter. The physical features vary only with the crossfrontal coordinate, y. In the linear stability analysis, the perturbation fields can be treated by separation of variables, and the parts in x and / are harmonic fields. The perturbed state is represented for each layer by four fields: Uj(x,y,z,t) = Uj + sUj{y)e«kx-"> + 0{f) Vj{x,y,z,t) = evftye***-* + 0(f) pj(x,y,z,t) = Pj(z) + £pj(y)e,(*x""> + Otf) hj(x,y,t) = h+y)* ( - i y etijiyyo*-* + 0(e>) where a: complex frequency (the real part is the physical frequency of a wave traveling in the ^-direction and the imaginary part is the temporal growth rate) k: wave number (assumed real) s. small parameter for the ratio of the scales of the initial disturbance and the basic state. 8 The variables Uj(y), Vj(y), Pj(y) are the first-order terms of the perturbation fields of alongfrontal velocity, crossfrontal velocity and pressure, respectively. The variable /7(y) is the first-order term of the perturbation in the thickness fields. It represents an anomalous elevation of the pycnocline. The governing equations in each layer are the nonlinear, vertically-averaged, inviscid shallow water equations without friction: du, _ du, „ du, -1 dp, —- + u —— + v ( .—-- iv, = -(1.4a) dt 3 dx 3 dy 3 p dx dv, „ Svj „ dv, -1 dp, —- + u,—- + v,.—- + fu, = -(1.4b) dt 1 dx ] dy 1 p dy dhj „ dhj m dhj ~ (1.4c) dt 1 dx 3 dy 3 (du. dv, J + 1 \dx dy J Eqns. 1.4a and 1.4b are the along- and crossridge momentum equations, respectively. Equation 1.4c is the continuity equation and h is same as hj in regions I and ///, but becomes the deformed layer thickness inside the frontal region (If). The separated fields are substituted into the system of governing equations, and the linearization is performed by removing zeroth order (8°) terms, that correspond to the initial balance, (1.1), and the products of perturbation elements that correspond to second order (s 2) and higher terms. Only first-order (e1) terms remain. Replacing the x- and ^ -derivative operators by the complex factors ik and /'a, the momentum equations become i{kU-o)u -Jv =-—p p / (^.-a)v y + / M y .=-I^- . p dy Let (Oj ={kUj -o)jf, (ratio of the Doppler-shifted frequency over the Coriolis parameter). The solutions for the velocity fields in terms of the pressure field are (15a) J (l-a?j)pf , (1.5b) v -to-*?* The subscript,y after the layer index represents differentiation with respect to y. The perturbed velocity fields can be determined only if the system does not oscillate at the inertial frequency. The continuity equation appears differently, according to the region. Outside of the frontal region (regions I and III respectively), the rigid lid approximation states that In the frontal region (IT), dhi dh. —L=o, thus —1=0-dt dx h\=H(l-y/L)-£T](x,y,t), h2 = h2(y)+ et](x,y,t). where et] is the vertical displacement of the pycnocline. The pressure of a water column is continuous at the pycnocline. The pressure under the pycno-cline is then (1.6) P2(y,z) + ep2 = Pl(y,H(l-y/ L)) - pgsrj+p2g(z - H{\ - y/L) + STJ). Removing the initial balance yields (1.7) Pg' The substantial derivative of the layers' thickness in region II yields: dt dx dt dx) pg' and dh2 dt dh dx drj , TT drj dx ) dt r io)2fs / x = +io)2fet] = +—J-r(p2 - a ) -pg In regions I and III, (1.4c) is reduced to vh =—h(iku+v ) j j,y y)ty} or, by substituting (1.5), (1.8) ; ~ »jPj^-j^°&hj(y» = -*>J(*2PJ - PJJ-In region II (the front), for the upper layer, (1.4c) becomes _ £ V i _ icoj?]^ -H(\-ylL){ikux + v j 10 or, using (1.5) and (1.7), 0-9) V H ( l ~ m ] )ip> - « ) = < » - y ' L * * * - " J + P i * ~ k L p ' " ° ' and for the lower layer, (1.4c) becomes viKy(<y)+iQ)2fTl=:-h2(y)(iku2 + v2y) or, using (1.5) and (1.7), (1.10) £ ( l - a>\){p2 - A ) = f P2,y ~ ^ k , 0 0 " ^ ( ^ A - Pl,yy\ The crossridge direction is made nondimensional by dividing by the frontal length: a i a y-> Ly, > , dy L dy thus U<-(l-c>))PfL> J (\-<o))pfL where K=kL. Let h^y) = i\(Ly) and h2(y) = /i 2(Zy)-The equation for outside the frontal region, (1.8), becomes (1.11) K K°>J •P,-P J dy The equations for inside the frontal region, (1.9) and (1.10), become (1.12a) ^H1~ - A ) = (1" .V)(* 2 A - A,J + P\,y ~ ~ A . V < W 2 A ~ Pj,y Let 11 A = £ £ . ( 1 - 0 ? ) , £=LL(i-ai)t ^A-K/^, girl g t i then (1.12a) becomes and (1.12b) becomes (1 -y)(fC2px -pXyy) + pXy + AxPx = Pxp2, 6), / / A = A A -These can be rearranged as (1.13a) (i-y) ( d2 „A d dy 1 Pi = -P\P%> (1.13b) H ( ^ dy2 -K2 H dy H Pi =fi'2Pi-This system of ordinary differential equations will use as boundary conditions the values of Pu Pi,y of region / aty=0, and of A>> Pt,y of region III at>>=0. 12 1.2. Triangular Ridge Consider the case of a triangular ridge (see fig. 1.2) with bottom located at H(\ - yxy) in region /, and at H(l + y2y) in regions II and III, where yx, y2, are nondimensionalised slopes i.e. slopes of the ridge divided by the pycnocline slope (H/L). In this section, we are going to solve for the pressure perturbation field in each layer and region by series expansions. x^r Region I Region II ^ v Region III Figure 1.2. Same as Fig. 1.1, except that layer thicknesses are defined as h\{y) = H(l- yxy) in region /, as h\{y) = H(\ + y^y) in region III, and, in region II, as h\(y) - H(l-y), h\{y) = H{\-*rY2)y for the upper and lower layers respectively. St] is the displacement of the perturbed pycnocline. Region /: In region /, the differential equation for the perturbation pressure, (1.11), becomes (1.14) CO, The problem can be simplified in two steps. First, change the independent variable y to z=l-y y; so (1.14) becomes Ar\p^-K2PxYy\pu+^-p^o. (1.15) Second, far away from the front, the perturbation pressure field must decay and become smooth. The balance inside the square brackets must vanish: 13 r2PU:-K2Pi<*0 as z->oo. Thus, the asymptotic behaviour of the solution is dominated by an exponential behaviour. Assuming that p, = Ajp(z, yx) exp(- Kzj yx), (1.15) becomes (1.16) 1-2K Yx '' J Pz 1 1 ^  CO p = 0. xj Making the substitution ax - \{\ - l/<y,) and the change of variable Zx = 2Kz/yx gives (1.17) ZxpZA + ( l - Zx )pZi - axp = 0, which is a confluent hypergeometric equation. Classical solutions of (1.17) are already known (see Erdelyi [1954]). The differential equation shows a singular point at Zx=0 (y = l/yx), which is not part of the region /. The global solution must behave in the far field such that it either decays smoothly or grows in some nonexponential manner: p(Zx)e -zxn - » 0 and p(Zx)e -ZJ2 ->0 as 00. This particular solution of (1.17) is known as the Tricomi solution: p(Zx) = W(ax,l;Zx). For a development that covers all the range beyond Z, = 2K/yx> it is numerically not practical to use a basis of the classical solutions. The solution is tediously exposed in Appendix A. It includes an asymptotic expansion from infinity, a sequence of intermediary matchings between series expansions (which can be replaced by a numerical solution of the ordinary differential equation) that extend the asymptotic solution to the point Z1 = 2K/yl. At this point, the solution provides two boundary conditions at7=0 for p,: (1.18a) A / ( 0 ) = ^ e - z ' / 2 ^ ( a 1 , l ; Z 1 ) | Z i = 2 J , / n and 14 ^(ax,\;Zx) i h-—(ax,l;Zx) dZt 'Z^lK/y, (1.18b) = AjKe -Z,/2 c9¥ V(axXZx)-2^-(ax,l;Zx) oZx Region ///: In a process similar to that of Region I, substitute 2F=l+7zy into the differential equation for Region III, which becomes tyiPi* ~ K 2 P i ] + f i P t , * ~ — P i = o. (1.19) 0)^ Assuming p2 = Amp(z, y2) exp(— Kzj y2), the ordinary differential equation (1.19) becomes ZPzz + K Pz n \ Ml) p = 0. The substitutions a2 =^(l + l/<y2) and Z2=2Kz/y2 give another confluent hypergeometric equation Z 2 P z 2 z 2 + ( l ~ Z2)pZi - a2p = 0. Its solution is found analogously to region I. It provides two boundary points at the endpoint of the front located at y = 1 (Z2 = 2K{\ + y2)/y2): (1.19a) / ? f (1) = Aine-z^(a2XZ2)\ _ 2K(ny2)/ri and p £ ( i ) = ^ dZ2 dy f V dZ, (1.19b) = -AmKe-z*12 f v F ( a 2 5 l ; Z 2 ) - 2 — ( a 2 , \ ; Z 2 ) oZ, J z2=2K{\+ri)/r2 'z2=2AT(l+r2)M Region II: For Region II, the Pi, Pi equations, (1.13a,b), form the coupled system (1.20a) and dl K 2 \ d dy2 ) dy P x = - P x P i 15 (1.20b) Let ( i + r 2 ) y dy2 -K2 dy K Pi + — 0 + r 2 ) CO, Pz = -p\P\-A* J2L2\-o?2 P l \ + y2 g'H l+y2 and (1.20b) becomes (1.21) 'd2 dy2 -K2 +• dy K co. 2 J Then substituting X2 = P2 + K/co2 and combining both equations gives dy2 -K2 + A, dy 2 (i-y) dy2 - K dy Pi =PAPI or, with D=d. 'dy [y{\-y)D*+(l-4y)D2-{2K2y(l-y) + 2 + X2 + y(XX -X2)}D2 + K*y(\ -y) + {{4y -l)K2 + X2 - XX }D (1.22) +{i + A2+Y(AX-X2)}K2+XXX2-fij32]px = 0. This leads to a fourth-order, ordinary differential equation for px with variable coefficients. It constitutes a boundary-value problem with the complex frequency as an eigenvalue, and two singular points at y=Q and y=l respectively. The physical boundary conditions assume continuity of A > A> > w i > w 2 > v i > v2 fields at the boundaries of region //: At JH>, A / (0) = A / ?(0) and /4(0) = /£ (0 ) . At y=l, Jpf(l) = pf(l) and /£(l)s=/>£,(l) . It will be assumed that regular singular conditions hold at the endpoints of region II. This implies that aty=0, (1.23) [D* - ( 2 + X 2 )D2 + (-K2 +X2-XX)D+(\ + X2 )K2 + XXX2 - pxp2 ]px = 0 and aty=0, (1.24) [ - 3 Z ) 3 - (2 + ^, )D2 + (3K2 + X2 - XX )D + ( l + XX )K2 + XXX2 - pxp2 ]p2 = 0. 16 1.3 The dispersion relation The Frobenius' method is used to find the series solutions for the frontal system for P l by expansion about y=0. Let (1.24) Then .n+r n=0 DPi(y) = lZan(n + r)y"+>-\ = 5 » + rXn + r - \)y""-\ ^ViW = + + r -!)(« + r - 2)(« + r - 3)y n=0 We can shift powers of>> to the same level by handling the indices. The various parts of (1.22) become y{\ - y)D4Pl + (1 - 4y)D'Pl = 5 > „ + 3 ( « + r + 3)(n + r + 2)(n + r + \)2yn+r n=-3 - ]T an+2 (n + r + 3)(« + r + 2)(n + r + + r)y"*r', » = - 2 - [2^1 - >0+2+v+^(i - y)p2p>=2*:22>>+'X*+^ - D / + r - (2K2+X-X2)Yjan+x{n + r +1)(« + r ) / * ' - (2 + 22) £ > w + 2 ( « + r +1)(« + r + 2)/+ r, n=-2 [(4j-l)K 2 -Al+A2]DPl = 4 ^ 2 + r)y"+r -(K2+ X-X2)^an+l(n + r + l)y"+r [y(l-y)KU ( ) * V +r] A= -K4^ a^'+T^ a ^ r «=2 M=0 +K2{K2U-A2)Yja„_ly"-where T = ( l + / l 2 ) ^ 2 + XXX2 - pxfi2 '2-The overall equation (1.22) becomes 17 £ an+i(n + r + 3)(« + r + 2)(n + r +1) 2f+r » = - 3 - X a „ + 2 [ ( « + r + 3)(« + r) + 2 + A 2 ](« + /• + 2X« + r + l)y"+r n=-2 - Z fl^i [(2» + 2r + I K 2 + (A, - A 2)(« + r +1)](« + r + + E a„ [2K2 ( » + + ' +1) + r Jy" + r n=l n=2 At each power level, the coefficients must cancel. (1.25) n = -3: a0r{r - l)(r - 2)2 = 0. In order to keep r meaningful, a0 must be arbitrary. Thus r=0, 1 or 2. There are only three roots of (125), then there must be only three fundamental solutions that have a regular singular behaviour about y=0. As the roots of (1.25) differ by an integer, their effect is to shift the power in the solution series. But it is sufficient to consider r=0, and to build fundamental solutions by factoring the arbitrary coefficients: n = -2: al(r + l ) r ( r -1 ) 2 -a 0 r{r - l)[(r + l)(r-2) + 2 + X2] = 0. For r=0, a, remains arbitrary. n = - 1 : a2(r + 2)(r + l ) r 2 = ^ ( r + l)[(r - l ) ( r + 2) + 2 + A 2 ] (1.26) + tf0>'[(2>* ~ ! ) ^ 2 + r U i " *i)] The case r=0 is again a trivial solution of (1.26), which leaves a2 arbitrary. For n>2, let k=n+3 (k>5) and use r=0 for the general indicial equation: k(k-3) + 2 + A2 (2k-5)K2 - v l 2 ) ( * - 2 ) " ^  * 0 T ~ 2 ) + a*"2 * ( * - l X * - 2 ) a t _ 3 [2£ 2 ( k - 2)(Jt - 3) + r] + ak_4K2[K2 + X, - A 2 ] - a A _ 5 £ 4 ( 1 2 ? ) A ( * - l ) ( * _ 2 ) 2 ' The coefficients a0, a,, a2 are arbitrary, so the set of fundamental solutions can be reduced to a linearly independent basis of three eigenfunctions, say: (1.28) 0 O O O = 1 - — y3 - ^ ^ J J - 1 _ J 2i / +... 18 with a0-l, a{=0, a2=0; KA+XX-X2 3 ( e+Ajf^+A. -Aj - fr + r ) v , . ~ y ,O0 = .y g y ^ y +••• with a0=0, <3,=1, a2=0; (1.30) 0 2 OO=v 2 — y 24 with a0=0, a,=0, a 2 - l . The missing eigenfunction that would close the basis of fundamental solutions is the one that contains the singular asymptotic behaviour about the regular singular point, y=0. It is rejected in order to respect the regular singular condition (the fourth derivative of px stays finite as y -> 0+). The global solution expanded about y=0 is (1.31) Pi =«o©oCv) + «i©iCv) + « 2 © 2 0 ; ) -The boundary conditions at y=0 imply 2K (1.32) a0=d(G) = Atx¥(alX—') and n (1. 3 3) = KA, OK OK ¥ ( a „ l ; — ) - 2 n « „ l ; - ) Yi Yi . where al, l,2K/yx) is solution of the transformed differential equation in the region I at the boundary, and W is the derivative with respect to the element after the semicolon, i.e., 2K/yx in this case. The coefficient az can not be determined on this side. Therefore, the global solution can be split into two solutions: one for which the slope and the magnitude at the origin depend upon conditions in the region J, and that has no inflection at the origin; the other solution behaves freely, without magnitude, nor slope, but curved at the origin. For p2, we get a similar fourth-order ordinary differential equation: [y(\-y)DA +(3-4.y)D 3-{2K 2y(l-y)+2 + X2 +y(Xx-X2)}D>2 +{(4y-3)K2 + X2-Xx}D+K4y(\-y)+ (1.34) K2{\ + X2 +y(Xx -X2)} + XXX2 -ft&h = 0. 19 Let (1.35) »=o The lowest power level in (1.34) is provided by the two first terms: [y(l-y)D4 +(3-430Z>3]/>2 = 6 0 ^ - l X ^ - 2 ) 2 / " 2 +0(/- '). So s= 0, 1, or 2. Similarly with px, after rejection of the singular solution, the series solutions can be treated as analytic and the three first constants {bx, b2, b3} remain arbitrary. Instead of creating another general indicial equation with five terms for the bn% we can use the fact that p , and p2 are coupled in order to shorten the algebra. Then, with the substitution of (1.24) and (1.35) in (1.20b), we get X k + 2 ( « + 2)(« +1)-aH+l(n + 1)2 - (A, + K2)an}y" + K2^an_xyn (1.44) - A l V ' Then, it follows from the regular singular condition , _ (A,x + K2 )a0 + ax — 2a2 and the first general indicial equation is a =a n-l (Xx+K2)an_2-K2a„_3-J3xbn.2 (1.45) " n n(n-\) With the substitution of (1.24) and (1.35) in (1.21), we get »=0 n=l n=0 Then *i =A2bQ-02ao and the second general indicial equation is (1.46) °« M 2 All 6„'s can be determined with respect to a0, a, and a 2. Setting a2 = 0 and prescribing a 0 and «j by the boundary conditions on the left side of the front, we get the parts of px and p2 20 that depend on the arbitrary magnitude of px in region I (Aj). Setting a0 = 0, ax =0, and leaving a2 arbitrary, we get the other parts of px and p2. Thus, A 0 0 = AjPx (a^c^, 0, K2,Xx,X2,px,p2\y)+a2Px (0,0,1, K2,XX, X2 ,fix ,fl2 ;y) p2 (y) = A{P2 (a0, ax, 0, K2, Xx, X2 ,px ,J32 ;y)+a2P2 (0,0,1, K2, Xx, X2 ,J3X ,/32;;;), where a0 = p((0)/Aj and ax = p^ify/A,. The radius of convergence of the series (1.24) is #c = l im—-=l im 7i->oo /7 , n-»oo n-l n + 0(n2) = 1. provided the nondimensional wavenumber K and the Xx, fix parameters stay bounded. And series (1.35) converges for any y if K and the parametersX2, ft2 stay bounded. Together, the series solutions do not converge at>*=0. We need another point from which the expansions can describe the neighborhood of y=l. Use the boundary point with region III and let z=\-y. Thus (1.20b) and (1.21) become d-4 (1.47) (1.48) (d2 t) d , —--K2 X2 Vdz J dz dz2 -K2 d , Pi - ~PzPi, Pi^-PiPt-The set of equations (1.47,1.48) with the boundary conditions from the region III stay similar to the original set, thus series solutions can be found similarly by swapping the indices: p2(z) = A'mPx(a'oXAK2,X2,Xx,02,^ px(z)=A\IIP2(al),a{,0X,X2,Xx,p2,px,z) + ^ where a> = pf ( 1 ) / ^ , a{ =-p™(\)/Am and 2Kh \ A'm - A, e 2 The solutions for px and p2 do not converge at z=l (jK)), but both sets can be connected at y=z=y2: ^ p 2 ( a 0 , « 1 , o , r ^ A 1 ^ 2 , A , A ; j ) + « 2 ^ ( o , o , i , A : 2 ^ 1 ^ 2 , A , A ; > ' ) = (1.49) ^ W , < 0 , ^ , ^ , ^ 1 , A , A ; z ) + ^ / » ( 0 , 0 , l , A : 2 , ^ 2 , A I > A f A ; z X 21 (1.50) ^ P 2 ( ^ , a I ' , 0 , ^ 2 f A 2 , A 1 , A . A ^ ) + ^2 ( 0 . 0 , l , ^ 2 , A 2 , A 1 , A . A ^ ) -The arbitrary constants of integration are A{, Am, a2 and a 2 . Two other equations are needed for the matching. They are provided by differentiation of the solutions evaluated aty=z=V2. AIP2XaQ,al,O,K\Al,A2,j3l,02-,y)+a^ (1.51) -^^, f l l ^o , ^^^ 2 ,A I f A,A^)-^Woa / :^A2^l .A .A;^ (1.52) -A\nP[(a^a[AK\X2^x,P2^ Primes denote differentiation with respect to y and z, according to the origins from which functions are expanded. The conditions that solve for the arbitrary constants form a linear where all the solutions are understood to be evaluated at y = z = j , with the parameters AT2, A,, A 2 , $ , /?2. The barred symbols refer to the swap of the indiced parameters Xx, X2, /?,, p2 in the independent solutions /?(•) and P2(-). The vector of nontrivial arbitrary constants lies in the nullspace of the linear system only for some specific values of the complex frequency, a, (or one of its avatars, such as the nondimensional, Doppler-shifted frequency, <a„ which forms a nonlinear eigenvalue). The condition for which a forms the eigenvalue of the fourth-order differential equation is that the determinant of the above system is zero. This condition forms the dispersion relation of the system. In the specific cases of fi2 = 0 or A* = 0 (the perturbation oscillates with an inertial period in the referential frame of one layer), the model loses its validity since the transformation from the perturbed velocity fields into the perturbed pressure field cannot be inverted. system: (1.53) i?(flb,fl„0) P,(0,0,1) -P2(a'0,a[,Q) -P2(0,0,1) P^a^O) P,'(0,0,1) PM,a[,0) P^O.0,1) P2(a0,ai,0) P2(0,0,1) -Px(a'0,a[,Qi) -/} (0,0,1) PXa0talt0) P2'(0,0,1) PXa'0,al,0) J>'(0,0,1) _ 22 14 Initial parameters The perturbation of the frontal system is described by a fourth-order differential equation with four boundary conditions; therefore the eigenvalue problem has four degrees of freedom. The set of parameters defined in this chapter must be constructed upon these degrees of freedom and will be used as original inputs. The first input parameter to consider is the nondimensional wavenumber K of the distur-bance. Considering that the system is frictionless, we have the freedom to choose the frame of inertia to match the centre of motion between the water masses. ,U,+U2 a Let co=k~ — 2/ / be the eigenvalue of the system. The next parameters will not require further knowledge of the momenta other than the momentum difference or nondimensional Doppler shift. From the eigenvalue, cox and co2 are readily defined: K{UX-U2) . K(UX-U2) co = CO+—1-1 IL and ay = oo —! £i . 2fL 2 2fL But Ux-U2 = AU = {AU)2 fL ~ fL~ g'H ' Let m_f2L*_ g'H g'H (At/) 2 " This new input parameter is analogous to a Richardson number and is a measure of the ratio of the available potential energy over the relative kinetic energy. Its reciprocal can be also seen as the squared ratio of the baroclinic Rossby radius over the frontal length, i.e. as a squared Burger number. It is the parameter F in the work of Gawarkiewicz [1991]. An alternative selection of input parameter would have been a Rossby number equal to the nondimensional Doppler shift, Ro = K/(2Ri), which was preferred by Kotschin [1932], Orlanski [1968] and Flagg and Beardsley [1978] instead of the nondimensional wavenumber. 23 Thus cox = co + K/(2Ri) and ct)2 = co - K/(2Ri) Four other parameters result from the above inputs and the eigenvalue. From the boundary conditions: " i 2 a. 1-1 CD 1 + — CO '2 J From the pycnocline region: px = Ri(\-co]l CO, The straight slope under the pycnocline calls for a third input, y2, and two other parameters: J32 = Ri A-co' 1+ y2 K co. In his treatment, Gawarkiewicz [1991] assumes perturbations close to a geostrophic balance, such that |of |, \a\\«1. From that approximation, fix and (1 + y2)fi2 a r e equivalent to his Richardson's parameter (g'H/(AU)2). So, f3x and fi2 are denoted modified Richardson's numbers. The last input to define is Yu and comes with the left boundary conditions. It was implicitly set to zero in the shelf-break model of Gawarkiewicz [1991]. The set of inputs, {yx,y2,K, Ri}, is made from real-valued parameters. Thus for any solu-tion of the dispersion relation, the complex conjugate of the eigenvalue is also solution. This implies the system is either unstable or marginally stable (no growth rate, no decay rate) for a set of original inputs, but can not be stable (in the sense of allowing only a decaying mode). 24 CHAPTER 2 NUMERICAL SOLUTIONS OF THE DISPERSION EQUATION The implicit dispersion relation for the eigenvalue is formed by the determinantal equa-tion of the system (1.53). It is solved numerically with a C program performing three distinct tasks: a computation for some point of the determinant on the left-hand side of (153), a search of its roots over some domain of the complex plane, and a decision if an apparent root is acceptable. At the lowest level, it will pick information about the configuration of the model topogra-phy and "Richardson's" number at the pycnocline), set files to record data, and run the search of eigenvalues for a spectrum of wavelengths with trial values (iterates) of the complex frequency picked from a specified region of the complex space. The method of search is by gridding a domain in the complex space and by identifying signs of the real and imaginary parts of the determinant equation over all the corners of every rectangular gridcell. When both signs change, contours of zero value of the real and imaginary parts (zero-lines) cross through the gridcell; the cell is checked with a bidimensional version of a bisection scheme. The determinantal equation of the system (1.53) is an implicit function for eigenvalues. For iterates of complex frequency, the left-hand side is not generally zero and can play the role of an objective function that will be hereafter named dispersion function. It is clear the dispersion function is not the dispersion relation which only eigenvalues solve. It is embedded in the C program through a function of complex type named dispersion. The dispersion function is a function of one complex variable and constant parameters. So it is holomorphic over the scanned domain. This is not always true for all points. There are two known points that constitute poles of the function {o) = (±Ro, 0)) and four other points for which the model does not hold (for /5,=0 and /52=0, the perturbed momentum fields can not be inferred from the perturbed pressure field and its derivative). For the analytical domain, we use 25 the Cauchy-Riemann theorem. Let w(x + iy) = u(x,y) + iv(x,y), where u(x,y) = Re(w(x + iy))t v(x,y) = Im(w(x + iy)) and w is analytic atx + iy then du _ dv du _ <?V Thus the contour lines of the real part are orthogonal to contour lines of the imaginary part at every analytical point. The parameters being real, the imaginary part of the dispersion function is odd with respect to the imaginary part of the iterate. Thus the real axis is a zero-line for the imaginary part of the function, and the zero-lines of the real part cross the real axis at a right angle at the points of marginal stability. Similarly, at points representing an unstable eigen-value, the zero-lines meet orthogonally and cannot be tangent. The size of the initial gridcell is assumed small enough so the zero-lines for the real and imaginary parts have very small curvature (sufficient resolution assumption). Therefore, neither can a contour line get into some cell and get out by the same edge, nor can it be confined inside one gridcell. A root involves the two zero-lines crossing through the initial gridcell, but the converse is not necessarily true: the root can be confined inside a near cell. The false root is detected by gridcell subdivision until the resolution sets the zero-lines apart. For the same assumption, a gridcell cannot have more than one root. Unfortunately, this assumption does not allow the routine to find closely paired modes at critical values of some parameters. Sufficient resolution does show locally orthogonal zero-lines about the roots, because numerical inaccuracies negligibly effect the analycity. But it has been observed untrue for extreme numerical resolutions: contours get smeared by noise. Nonetheless, the root is accepted. Moreover, there can be points for which some eigensolutions become either unbounded or the basis of fundamental solutions loses its linear independence before matching. In the latter case, the subdeterminants between the entries of free and boundary-tied solutions from the same side of the front are numerically negligible with respect to those entries. This can be caused by exploding solutions. It becomes a stiff problem at those points and spurious roots appear due to computational noise. In all those cases, the routine dispersion is designed to detect the problem and return a flag to signal it. 26 The determinant of (1.53) is computed as «„ W 1 3 « 1 4 ™21 W22 m23 /w24 W31 ™32 « 3 3 W34 w4 2 W43 » » 4 4 m12 W33 ™34 + ™32 ™23 W24 /w12 m23 w 2 4 w2 1 m22 m43 «44 mn W43 W44 W 4 , W42 W33 ™34 m2l m22 m13 + m42 m,, mu + W31 >»32 mi3 m14 m3l m32 m43 w4 4 m2l m22 W33 »»34 »»41 W42 m23 «24 Let ^ 1 £/3 = t/5 = mn mn > A = W33 mj4 , u2 = W32 . A = W23 ™24 m22 ™43 W44 TO,, >"l2 m43 ™44 mn ml2 , D 3 = ™23 »*24 . ^4 = w21 2^2 . A * « 1 3 * I 4 m4l mn W33 W34 W31 »*32 »»43 W44 m41 m42 > A = W13 m14 " » J I ^2 . A = 1^3 « I 4 m21 m22 m33 >»34 « 4 1 W42 W23 ™24 If, using £• as a criterion of numerical smallness, max« m a x ^ ^ l ^ i 1 ^ } ' max{|w,,m 2 |, |»i l lm n |}' max{|w11w42|,|7«417«I2|}' M N M max{|m„»% 2 m i i« a |} ' m a x ^ ^ l ^ i / n ^ l } ' max{|/^,w42|,|?M41m32|} i.e., it is relatively negligible for small 8, then the basis of numerical eigensolutions on the right-hand side of the matrix (free and bounded solutions expanded from y=0) loses its linear independance. That indicates a stiff problem. The same criterion is applied on the left-hand side. Calculations at this level are done in double precision and E is set to 10"9. The determinant is the dot product between the }and {Dy} vectors. The neighborhood of a root may consist of terms that could cause overflow after conversion in single precision. The determinant is renormalized in its real and imaginary parts by respective greatest element before last addition. f det = max lRe(t/,D,) max JIm(C/,D,) v>{l,2,3,4,5,6}l 1 J \ y={l,2,3,4,5.6}l ^ •'I The subroutines are described in detail in the appendix B. 27 CHAPTER 3 RESULTS 3.1. Flat cases 3.1.1 Low Richardson number The program was first tested for the case of the Kotschin model investigated by Orlanski [1968], that is a flat topography with a weak Richardson number (Ri=3), applicable in meteor-ology. Orlanski [1968] plotted growth rates against Rossby numbers. He found two pure imaginary eigenvalues for small wavenumbers (K< 0.1) that coalesce at Ro=0.1 to give birth to a pair of symmetric unstable modes (same growth rate, opposite frequencies). According to his published plot, which is a continuous curve, the growth rates peak at Ro=0.25 with a value near 0.107, then fall to zero tangentially to the axis when Ro=0.72. For higher Rossby numbers, another instability appears and peaks at 0.08 at about i?o=1.12. His curve appears to be an interpolation based on 24 points. Following the points from his appendix C, his curve cannot be very accurate: most points describe the first peak and the highest growth rate is at Ro=A, with a value of 0.133, while the second peak is not as smooth as shown on the curve. The double instability at small wave numbers is not present in the observed points for 7?/-3, but is extrapolated from Ri=2.2. Our results (Fig. 3.1) show a very good agreement for the shape at small Rossby numbers and Orlanski's points (not shown) fit well on the first peak. However, here a cutoff appears before Ro^O.7, with everything beyond stable, whereas Orlanski [1968] only found local stability in this region of Ro and further instability at higher Ro. With the transformation r= co/Ro, we improve the description of unstable features (see fig. 3.2). The double instability starts from nonzero points: the upper branch starts from x=(0,l); the lower branch starts from T=(0,0.4) as computed by Orlanski [1968] for very long waves. The upper branch is part of a surface denoted "R" by Orlanski [1968] that continues a Rayleigh-type instability (shear waves) deformed by the Earth rotation and the frontal geometry. The lower branch is part of his E surface. He characterized the upper branch as of barotropic 28 instability because that mode feeds off the mean kinetic energy, while the lower branch was characterized as of baroclinic instability because it feeds off the mean potential energy. F l a t t o p o g r a p h y : ( O r l a n s k i c a s e ) R i = 3 F = . 5 7 7 3 5 2 . 0 1 . 0 0.0 R o s s b y n u m b e r r a n g e ( R o ) : P h a s e f r e q u e n c y , „ , ! ! ! ! : ! ! " 3 a o E O . 1 0 G r o w t h r a t e 'I I it I y 0 1 2 3 4 5 6 7 8 9 1 0 1 1 nondimensional wavenumber K Figure 3.1. Flat bottom and Ri=3. The upper graph represents a nondimensional Doppler-shifted frequency against wave number. All modes are contained inside the range (-Ro, Ro) shown by dashed lines. Due to symmetry: only the positive frequencies are shown. The lower graph is the nondimensional growth rate. At Ro=Q.l, the double instability branches join to form a new unstable mode with nonzero real component and the imaginary part decreasing with a linear slope of 0.9 relatively to the Rossby number scale until Ro=QA. This slope gives a parabolic shape to the growth rate peak. The new unstable mode is part of Orlanski's [1968] B surface that is a mixed type (baroclinic-barotropic) instability, but in fact is a baroclinic instability at this Richardson number. After the cutoff point, he predicted a Helmholtz-type instability (H surface) deformed by the Earth's rotation and frontal geometry. This kind of instability did not appear in this present study. Orlanski [1968] mentioned that his technique converged poorly for Ri>2 and Ro < 1, which may explain the discrepancy. 29 F l a t t o p o g r a p h y : ( O r l a n s k i c a s e ) R j = 3 F = . 5 7 7 3 5 1 ' ° r P h a s e f r e q u e n c y / R o s s b y n u m b e r o 0.9 0 .8 0 .7 0 .6 0 .5 0 .4 0.3 0.2 0.1 0 .0 stable stable unstable • stable 3table 0 .0 0.1 0 .2 0 .3 0 .4 0 .5 0.6 0.7 O.f _ G r o w t h r a t e / R o s s b y n u m b e r o 3 0 .8 0 .6 0 .4 0 .2 0.9 1.0 1.1 1.2 1.3 1.4 1.5 R o R V 0.0 0 .4 0 .5 0 .6 0.7 0.1 0 .2 0 .3 R o s s b y n u m b e r : R o Figure 3.2. Flat bottom and Ri=3. The Orlanski's eigenvalue representation with Ro as a state variable. Due to symmetry, only the positive side of the real part is shown. Note the horizontal axes have different scales. The symmetry of the equations relative to the eigenvalue in the flat bottom model case allows the implicit dispersion function to be even. That symmetry is carried by the real component of the eigenvalue, for which only the positive part is shown in figure 3.2. In the range Ro<0.l, two unstable modes of different growth rates share same frequency with Re(<y) = 0. Stable modes start horizontally from t=<±0.59,0), (±0.79,0), and (±0.873,0). Other stable modes might exist closer to the Rossby barriers (co^iRo), but the vicinity of a pole makes them beyond the arbitrary safety range. The stable modes "feel" the end of the double instability as evidenced by a slight curvature toward the Rossby barriers. The continuation of B instability by stable branches forces a compression of the stable modes against the Rossby barriers. Hereinafter, the stable branches starting from the origin and staying close to the 30 Rossby barriers are denoted marginal stable branches. Following the growth rate peak at i?o=0.4, the real components of the unstable modes bend to become horizontal at about ±0.4 at i?o=0.6. That is approximately the cutoff point of the unstable range (A=3.8). Each unstable branch is followed by two stable branches forming a rightward parabola. While the outward branches compress the marginal stable modes, which tend to become parallel to the Rossby barriers for high wavenumbers, the inward branches meet on the axis and cross without apparent interaction at Ro=0.95. The computational range of the program is quite limited beyond i?o=1.15 (AMO), so we cannot observe their behaviour for large Ro. 3.1.2 High Richardson number Although the model for an oceanographic context is similar to an atmospheric context, the Richardson number is generally much higher. In agreement with parameters from Flagg and Beardsley [1978] and Gawarkiewicz [1990], it is set to 57. Extrapolation from Orlanski's [1968] instability surfaces suggests that the H surface borders the B surface as Ro->\, that is £=114, which is out of the computational range. It is not clear if the line separating the R, B and E sur-faces stays near Ro=0.1 or is confined to smaller values. For a spectrum extending to about A M 1.5, the transition is expected to be observable. Our findings (see fig. 3.3) for a range of wavenumbers between 0.2 and 10, show a different perspective. At high Ri, the eigensolutions for frequencies near Ro, behave similarly to each other and there are more eigenvalues. The density of marginal stable branches increases dramatically. The regions 0.75ito<|a>|</to are not graphed for clarity. The pattern of unstable modes has eight distinct branches. Their separation in the frequency component is largely proportional to the wave number K. Two branches (a ,e) are purely imaginary. The less unstable branch (e) has its maximum about £=1.75 and decreases to zero at AT=2.75. It gives birth to two stable branches tangent to the axis. The other branch with purely imaginary roots (a) reaches its maximum near A=7.4, © =0.02. Al l other unstable branches have greater growth rates than these purely imaginary branches. The outermost (largest phase frequency) branches 31 (d) have growth rates comparable to the major central branch (a) until K=6, where their growth rates keep increasing while the central mode growth rate starts to decay. The innermost branches (b) have higher growth rates than the major central branch and reach their maximum at about AT=9. The growth rates of the middle branches (c) are similar to the innermost branches (b) for moderate wave numbers, but they give the most unstable waves, at about AT=10, with a growth rate of 0.029. The outermost locus ends with negligible curvature, which hints at much less stable waves beyond AT=12. There are no stable points found outside the region -i?o<Re(o))</?o, hereinafter denoted Rossby range. Rossby number range (Ro) : •-Phase frequency of stable modes 2 3 4 5 6 Stable frequency/ Rossby number Phase f requency b a 0 1 2 0 0 3 f Growth rate • t r < a " ,,e K 0 1 2 3 4 5 6 7 8 9 nondimensional wavenumber K Figure 3.3. Flat bottom configuration with a Richardson number set at 57. On the left panel: the observed stable points (up) and their t-transformation (down). On the right panel: the observed unstable points with frequency (up) and growth rate (down). Due to symmetry, only positive real parts are shown. Reusing the transformation x=a>/Ro improves the description for small wave numbers (see fig. 3.4). Although it does not provide clear picture for what happens at tiny wave numbers, the branches seem to have finite ratios of growth rate over Rossby number as AT->0. The connection of the purely imaginary branches with the R and E surfaces cannot be established. The R mode should have its intercept at x=l, while the E mode should have it about xi=0.29 32 according to Orlanski's [1968] fitting formula; although the intercepts of the imaginary branches are not clearly defined, they do not follow the expected values. In addition, they do not join in a B surface of only one unstable mode starting at their junction. Instead, there is more than one branch with complex roots spanning over all the visited wavenumbers. While the inner branches (b, c) have similar growth rates for moderate wave numbers, the outermost branches (d) have similar growth rates as the central branch (a), but with opposite curvature. The central branch (a) tends to a cutoff at about A M 1.5 or i?o=0.1. Phase frequency/Rossby number 0.8 h o or 3 0.6 O ® 0.4 h 0.2 h 0.0 d H M 1 * " • * I • I I I t • I • I I I I I , | • l l t l l l l l l l I I » I I I I I M I I I I ! I i i i i c 11 • I I 1 I I 0 1 2 3 4 5 Growth rate/Rossby number 10 11 0.8 h 3 0.6 h ^ u  fo c a 0.4 0.2 h 0.0 .^.f 1 1"".'. '.-' n . . . m : . c , ' • . a b 0 1 2 3 4 5 6 7 6 9 10 11 nondimensional wavenumber K Figure 3.4. The t-representation at ^?/=57 for the flat bottom configuration. 33 3.2. Flagg and Beardsley shelf-break model The Flagg and Bearsdley model or false bottom model (hereinafter referred to as FB) consists of adding a straight sloping bottom beneath the front but the topography is flat outside this region. The slope y2 (equivalent to R in Flagg and Bearsdley [1978]) varies between 0.25 and 2.5. It breaks the symmetry of the phase frequency and creates stable modes in the negative region outside the Rossby range. For /2=0.25 (see fig. 3.5), two well-defined stable modes exist outside the Rossby range with opposite characteristics. The first stable mode starts from the origin with strong negative group velocity (in the reference of the momentum center of the layers) for small wavenumbers (A!"«l). The group velocity wanes as the wave number increases and, ultimately, becomes zero before the Rossby barrier. The second stable mode starts with positive group velocity and curvature from a finite intercept. It reaches zero group velocity at about A=3, then the branch straightens gradually and tends to join the Rossby barrier much beyond the observed range. Inside the Rossby range, the symmetry of stable branches is broken by two sigmoidal branches starting from the upper region to join the lower Rossby barrier. The marginal and central stable branches observed in the flat bottom configuration maintain a symmetric pattern. The main change is the beginning of the central branches is moved to about A=5.5 and the lowest marginal branches are moved slightly lower for the small wavenumbers range. The central stable branches are hereinafter labelled parabolic branches for their shape when unified in the observed range. Several distortions occur in the marginal branches where the sigmoidal branches cross. For the unstable modes, loci can be regrouped in six branches labeled A, B, C, D, E and B'. The growth rates are generally constrained beneath /to/3. The A, B, C, D branches dominate with nearly the same growth rate in the high wavenumbers range- Flutterings in the growth rate curves for greater wave numbers suggest slightly different growth rates for same frequency at same wavenumber, as in the branch C. The branch C is the only one converging asymptotically toward the axis for high wavenumbers and reaches a maximum at about £=8.5. The outer 34 branches, A and D, have the greatest growth rates beyond A=8, (greatest observed growth rate is about 0.025). The branch B seems to be an extension of the branch B' in the positive frequency region. The branch E, with eigenvalues nearly imaginary, reaches a maximum at about 0.005 within 3<K<4, then a cutoff happens about A=5.5, where it is followed by the two parabolic branches. The branch E is the branch e from the flat bottom configuration that becomes more unstable and extents to higher wavenumbers. Extrapolation of the real and imaginary parts of the branch D suggests that it may originate from the lower Rossby barrier, near AT=3. 0 04 r P h a s e f r e q u e n c y +R0 o.oo Doppler—shifted frequency -Ro -0.01 -0.02 -0.03 -0.04 0.03 f-0.02 001 f-, / ' . 5 5 M l • ' " 0.00 M>»S8J&3 43533* 2 2 5 2 2 2 2 2 2 2 E B' , , 1 I I 1 1 1 111 1M 1< ' C i 5 S t B 4 4 4 4 < , -Ro T = W / R O 0 1 2 io 0 0 3 r Growth rate i * t i i i i i t i t i , • i i i i t i * i i i i i i i i i , I • i i i / Ro A J 2 E 2 K K Figure 3.5. FB configuration with y = 0.25. Left panel: the frequency of stable modes outside the Rossby range (up) and the t-transform of stable branches inside the Rossby range (down). The dashed line is the lower Rossby barrier. Right upper panel: the real part of the unstable points. Dashed lines are the upper and lower Rossby barriers. Right lower panel: the imaginary part of the unstable points. For ^2=0.4 (fig. 3.6), the separation between the stable branches outside the Rossby range is enhanced by the stronger slope. The most unstable point was found at the end of the branch D with a growth rate of 0.0218. The dominant branches share roughly the same growth rate curves: discrepancies in the high wavenumbers range are mainly due to fluttering. The branch D keeps a fairly constant value of frequency at about -0.022: its apparition after K=5 forces the branch C to 35 converge toward the axis. The branch E maximum moves to AT=4.8, but the branch terminates about A>5.5. The point of the maximum in branch E can be seen as a crossing point between branch E and branch B. The dominant growth rate curves are negligibly affected by the increase of the slope, apart from a slight depression for the moderate wavenumbers. Inside the Rossby range, a new stable mode appears for strong negative values of t at small wave numbers. The central parts of the sigmoidal branches move to higher wavenumbers and the parabolic branches start about K=6.25. 0.00 to Doppler-shifted frequency - R o 0 0 4 r Phase frequency 0.03 |. 0.02 o.oi -o .o i -0.02 -0.03 + Ro » « 0.00 k » ' W " . , 3 » ! ! S 2 22 2 » J " < ' [ X . " 1 11 i n i i m n " ' ' . 3 E 2 2 2 2 2 2 2 B , ] i H i l l I 3 3 C ] J . 5 s s D -Ro = u / R o 10 0.03 Growth rate Ro . i t i s 14 <«*' 9 K 1 0 Figure 3.6 Same as fig. 3.5, but for a FB model with y2=OA0. For /2=0.63 (fig. 3.7), the most prominent feature outside the Rossby range is the apparition of another stable locus along the lower Rossby barrier as the first stable mode expands away. The major change for the unstable branches is vanishing of the branch D which allows branches E and C to reach more negative frequencies for moderate wavenumbers. The branch E reaches a maximum growth rate at AT=6, and its cutoff at K=H. Inside the Rossby range, the stable branch, that first appeared for y2=0.4 with negative frequencies only, extends to higher wavenumbers and is on the brink of positive frequency at very small wavenumber. The start of the parabolic branches suggests the cutoff of branch E occurs in 7<Af<7.25. The 36 distortions of the marginal branches are enhanced by the higher slope: most crossings between the upper marginal branches and the sigmoidal branches vanish. Instead, every sigmoidal branch is deflected rightward onto the former path of the marginal branch. The marginal branch is itself deflected upward to distort the next marginal branch and take its path. The result is a chain of distortions on the original path of the sigmoidal branches. In the lower region, the sigmoidal branches cross two marginal branches before being deflected rightward and inserted as a new marginal branch. Doppler-shifted frequency 0.00 U -0.05 0.04 0.03 |-- R 0 0-02 0.01 0.00 -0.01 -0.02 -0.03 \--0.04 Phase f r e q u e n c y ,--"' + Ro J * ^ 1 3 » « 4 4 4 4 * * * 1 4 ' i.t^f'oi* 4 44 4 4 4 4 4 , , \ . 1 1 » m m . m i l l 1 " " E - R o , „ 4«44*M4 « * * B mxt a »a' C ! 1— 1 2 T = W / R O 1 i : Growth ra te Ro IT A * i 8 A* *4 \ » Wftl 41* 1 1 * 1 V « * 4 4 I t , ,34>J 1** 1 E 1 9 K '° Figure 3.7 Same as fig. 3.5, but for a FB model with y2=0.63. For y2=\.Q (fig. 3.8), the new stable branch, observed for the case /2=0.63 outside the Rossby range, expands below the lower Rossby barrier. The branch E moves into the previous position of the branch C in the moderate wavenumber range but merges with the branch B before K=S. The branch C seems to have moved its root from the origin to some moderate wavenumber on the lower Rossby barrier. The main feature of the unstable branches is a general drop of about 0.005 in the growth rates beyond K=6. The parabolic branches inside the Rossby range start at A=8 but do not continue the unstable branch E. 37 Doppler-shifted frequency -Ro 0.04 0.03 |-0.02 0.01 0.00 -0.01 -0.02 -0.03 I--0.04 Phase f r e q u e n c y + Ro 44 4 « 3SD 3 S 333 -Ro 2 » S 2 2 C 3 4 5 T = ( j / R o : , i Growth ra te Ro 2 C 33* *4* 1 * 3 4 5 6 7 —1 \— 4 5 Figure 3.8 Same as fig. 3.5, but for a FB model with x2=1.0. For y2=\.6 (fig. 3.9), a new unstable branch, denoted F, appears with a positive frequency about two thirds of Ro and is much stabler than other unstable branches. Although this branch might have been masked by the upper Rossby barrier under smaller slopes, it might also be a set of unstable points induced by the large distortions on the upper marginal stable branches. The branch D reappears rooted in the Rossby barrier beyond A7=8, but it might be induced by the distortions in the lower marginal stable branches. Branch E appears to merge with branch B at about A>9. 38 Doppler-shifted frequency -Ro 0.04 R 0.03 0.02 0.01 0-00 < 4 4 • U ; : 1 4 " ! 4 4 4 - 4 4 4 4, * « * -O .01 - 0 . 0 2 - 0 . 0 3 - 0 . 0 4 - 0 . 0 5 - 0 . 0 6 - 0 . 0 7 Phase f r e q u e n c y + Ro F 111 1 " 1 4 B - R o T = W / R O 8 9 0 - 0 3 r Growth ra te Ro • 4 « 1 ? 3 * * ' J C 5 S D Figure 3.9 Same as fig. 3.5, but for a FB model with y2=l.6. For /2=2.5 (fig. 3.10), a fourth stable mode appears beneath the lower Rossby Barrier. Branches C and D have vanished. Branch F is enhanced and moved further away from the upper Rossby barrier. Its growth rate curve, although comparable to the other growth rates, suggests cutoffs at about K=4, and beyond AT=7. Branch A undergoes a transition near A>8. Branch E appears to coalesce with branch B at about A>9.5. The growth rates are generally constrained under the line Ro/S. The parabolic branches, starting at K=9.5, have almost moved out of the observed range of wavenumbers. Instead, the distorsions generate a new type of branch from the union of the least (smallest x) marginal branches and the sigmoidal branches. 39 Doppler—shifted frequency -Ro 0.O4 - 0.03 I-0.02 0.01 0.00 -0,01 -0.O2 -0 .03 -0.04 Phase frequency + Ro 2 22 -Ro T = O J / R O 5 6 7 i "—i 1 1 1 r— r 0 1 2 3 4 5 6 7 8 9 10 9 ° ' M f Growth rote / Ro • 3 F 5 6 7 Figure 3.10 Same as fig. 3.5, but for a FB model with ^ 2=2.5. In general, the real part of the unstable branches were constrained to where marginal branches are absent. As the marginal branches approach the center with an increasing slope, distortions change their nature and allow the unstable branch F in their region. The increasing slope causes a migration of the intercepts of the sigmoidal and parabolic branches to higher values on the iT-axis. With the exception of branch A, the unstable branches migrate toward negative frequencies: the branch C replaces the branch D; the branch E moves to the C place except for a part attached to the branch B; most of the branch B becomes purely imaginary instead of the branch E. Inside the Rossby range, the wavenumber where the parabolic branches start (Kc), varies monotonically with the slope (see figure 3.11). Curve fitting on the estimates yields: Kc = 4.62784 log(44.699 y2 + 6.36277). 40 io r 7 -6 -.4-5 - / F i t t e d c u r v e K c = 4.628 tag(44.7 7 2+6.38277) A -/' , 1 1 , 1 , , 1 1 , — 1 0.0 0 .5 1.0 1.5 2 . 0 2 .5 7 2 3 .0 Figure 3.11 Critical wave number for parabolic stable branches Kc vs y2. The stable modes outside the Rossby range are classified in two types by their behaviour in the long wave range. The first type covers modes with negative group velocity of norm l higher than YkT> ^ u t m e frequency is zero in the limit of long waves. The second type is for the backward propagating mode with finite nonzero frequency for very long waves. A third type l could be for modes with positive group velocities higher than YrT m& z e r o frequency at the limit of long waves, but it was not observed. Table 1 contains eye-extrapolated values of group velocities and intercepts for the fastest mode of the first kind and the mode of second kind at different values of y2 (graphic extrapolation made from points at small wave numbers). Intercepts of the first type modes are close to the origin, so the separation of the intercepts at A>0 is basically the norm of the intercepts for the second type mode. r2 Vg first type V g second type intercept separation .25 -.026 .022 .111 .40 -.037 .036 .165 .63 -.050 .056 .239 1.0 -.071 .077 .333 1.6 -.097 .093 .442 2.5 -.125 .096 .556 4.0 -.149 .097 .666 41 Table 1. Group velocities of the fastest mode of first type and of the second type mode interpolated for very long waves w.r.t. the slope. The intercepts' separation is the extrapolated value of the frequency for the second mode for long waves. The points seem to behave smoothly (see fig. 3.18, page 47). Using trigonometric and hyperbolic functions modeling with three parameters, the best fits are (r2+.1608) first type fastest mode: - . 2826 arctan second type mode: . 5991 arctan 1.475 (/2+.1137) intercept separation: .5580 arctan 1.021 (r 2 +1015) 1.655 42 3.3. Shelf-break model with foot at y-2 The change for this model with respect to the FB model is that the boundary conditions of an exponentially decaying perturbation pressure occur not at the front but at the foot of the submarine ridge and the dynamics between the front and the ridgefoot obeys a hypergeometric differential equation. In the following cases for various slopes, the stable branches beneath the lower Rossby barrier increase like interference patterns in optics where the slope, y2, is analogous to a lens curvature. The symmetric pattern of parabolic and marginal branches inside the Rossby range stands firmly, although the marginal branches are distorted by the skew stable modes. The most prominent change is the reshaping of the sigmoidal branches: their top part turns gradually rightward and there are smaller distortions on the marginal branches. A consequence is a reduction of skewness in the distortions relative to the FB model. Generally, the first skew branch on the right is nearly shaped as a parabola. The next skew branches are vertical in the central range, which make them nearly unobservable, and turn rightward in the region of the lower marginal branches. They migrate rightward as the slope increases. °'041" Phase frequency + R 0 i i u i i in i i i i C Figure 3.12 Same as fig. 3.5, but for a shelfbreak with y2=0.25, Ri=57 and foot at y=2. 43 For y2=-25 (fig. 3.12), a third stable branch appears beneath the lower Rossby barrier. The intercept of the highest-frequency stable branch is farther from the origin. The longer slope creates a new unstable branch, denoted G, which appears to be rooted on the lower Rossby barrier at about K=4. Branch G shares a strong similarity with branch D, which appears to be rooted on the lower Rossby barrier at about K=2.5. The growth rates of branch G are considerably smaller than those of the other unstable branches, but the difference vanishes for high wavenumbers. Branch E is well separated from branch B. Other unstable branches did not change much. For y2-A (see fig. 3.13), branch G has vanished. It is not clear whether branch B could meet branch E at about AT=4.5 as branch B's frequency changes sign. Coincidently, the major branches are slighthly more unstable between K=2 and K=6. The growth rates curves are still constrained beneath the line Ro/3. No other significant changes occur. o w r Phase frequency , - ' + R o - R o 0.03 0.02 0.01 0.00 -0.01 l--0 .02 -0 .03 i > ? \ I 1 ? 2 1 " ' , 1 i n i l l ' ' I i l i t i ' 0 I ! S S l l S i S I C x , -Ro T = C J / R O -0 .04 0 1 2 3 4 5 a M r Growth rate / Ro 3 . s ' 5 A 3 4 » ° I « ' 3 ' , S 3 3 5 , 1 * , ' i 2 —-1 —\ 1 ' 1 i 1 1 1 r— 0 1 2 3 4 5 6 7 8 9 Figure 3.13 Same as fig. 3.5, but for a shelfbreak with /2=OA0, Ri=57 and foot at y=2. For /2=.63 (see fig. 3.14), another stable branch appears along the lower Rossby barrier. Meanwhile, a new skew branch starts in the region of the upper marginal branches at small wavenumbers. At the exception of the branch A, the growth rate curves are depressed at about 44 K=6. Branch D is still observed at this slope, in contrast to the FB model case, and seems to be rooted at about K=5. Branch E coalesces with branch B at about A==6.5. The beginning of branch C shows a decreasing growth rate with K as it reaches a maximum negative frequency. Meanwhile, branch E looks perturbed at about K=3, with a kink in its growth rate: it may suggest a fork and a possibility would be that the unbserved part would connect to branch C. Branches A and B are not significantly changed. o.oo r-i.^--0.05 -0.10 -0.15 -0.20 -0.25 -0.30 -0.35 -0.40 -0.45 0 0 4 (" Phase frequency - R o 0.03 0.02 0.01 I-0.00 -0.01 |--0 .02 - 0 . 0 3 +Ro V 5 3 3 1 1 1 1 1 1 s 7 C I 2 2 A i 1 B 5 3 5 5 5 4 4 . D , -0 .04 0 1 2 3 ° ' 0 3 r Growth rate -Ro Ro 2 2 3 « 5 2*1*8 5 C 5 ' * 2 5 J 4 ' . 2 2 ? 2 5 4 E ! ! . / , , 3 6 6 , ' ? 1 » I 3 1 ' ~! 1 ~r-1 2 3 9 K 1 0 Figure 3.14 Same as fig. 3.5, but for a shelfbreak with y2=0.63, Ri=57 and foot at y=2. For x2=1.0 (see fig. 3.15), branch D vanished. Neither branch C, nor branch B changed with respect to the FB model. Branch A remains largely unchanged with respect to the previous slope, but its growth rate curve does not collapse at about A=6 as in the FB case, instead it follows largely the line Ro/4. Branch B remains largely similar to the FB case and coalesces with branch E at the same location. Branch E is finally broken by the kink observed at the previous slope; the beginning of the main branch is moved away in the negative frequencies domain and undergoes a hike of its growth rate. The minor locus, labelled E', terminates at K=4, but hints for a near-coalescence with branch B at about K=5. As the skew stable branches move 45 to higher wavenumber inside the Rossby range, their central part get rounded. The marginal branches appear more distorted at smaller wavenumbers by the skew branches. °- 0 4 f Phase frequency + R 0 o.oo -0.05 -0.15 -0.20 -0.25 -0.30 -0.35 -0.40 -0.45 -0.50 -Ro • 0.03 0.02 0.01 0.00 -0.01 -0 .02 -0 .03 9 -0 .04 i i 1 1 1 1 , 1 1 ' B 4 4 j E 1 1 « » 6 -Ro = w / R o 0 1 2 3 4 5 0 0 3 r Growth rate / Ro i 1—•—i 8 9 10 16*48 1 c ? 1 . 8 1 c E -1 !— S 9 Figure 3.15 Same as fig. 3.5, but for a shelfbreak with y2=l.O, Ri=57 and foot at y=2. For y'2=1.6 (see fig. 3.16), branches D and F from the FB model are absent. Branch B suffers two slight distortions relative to the FB model. Branch E is split near A=4 as in the previous slope, but the growth rates peak is absent. Branch E is disturbed again at about K=6 with a surge in the growth rate curve. The same occurs on a smaller extent to branch A at about K=6. 46 - R o 0.04 0.03 0.02 0.01 0.00 -0.01 -0.02 -0 .03 Phase frequency +Ro 13 3 3 3 3 J j j 3 3 3 3 3 3 3 3 , , .1 ' i t , 1 1 1 ' , 1 ' x - R o B 2 . 2 * E E 2 4 C 4 0.75 0 1 2 0 0 3 r Growth rate ' + R o 4 * C . 3 j 3 3 " 1 S 3 5 2 2 2 2 2 2 E 9 K 1 0 Figure 3.16 Same as fig. 3.5, but for a shelfbreak with y2=l.6, Ri=57 and foot at y=2. For y'2=2.5 (see fig. 3.17), branch F from the FB model is absent. Branch A is crossed by two nearby loci. The first locus, A', shares the same frequencies but differs by having a greater growth rate until they coalesce at £=3. The second locus, A", shares the same growth rate with A, but its positive frequency is greater until they nearly coalesce at £=3.5. Branch B suffers a strong frequency distortion, at its peak growth rate at £=5.5. It is noted that at B's end, points with the same frequency but slightly different growth rates appear; a lack of resolution may have concealed B's continuation. Only the part of branch E beyond its growth rate peak at £==6 is captured; its beginning point at £=6.25 provides the most unstable eigenvalue. Only two points of branch C remain. It is notable that the most unstable eigenvalues for these configuration parameters are produced by critical points at moderate wave numbers. 47 _-Ro 0.04 r 0.03 0.02 h 0.01 0.00 -0.01 -0.02 -0.03 -0.04 Phase frequency + R 0 Si,, H " J 3 l i t i 1 3 3 I ' ' 1 . . -Ro 1 M B T=6j/Ro —i 1 r~ 1 1 1 5 6 7 8 9 10 Growth rate Ro 1 . 3 > 3 A" ' 3 2 j 3 A 0 1 2 3 4 5 Figure 3.17 Same as fig. 3.5, but for a shelfbreak with /2-2.5, Ri=57 and foot at y=2. Vg first type Vg second type intercepts separation .25 -.072 .1025 .199 .40 -.102 .144 .283 .63 -.144 .175 .283 1.0 -.192 .213 .498 1.6 -.240 .238 .613 2.5 -.30 .225 .715 4.0 -.35 .200 .800 Table 2. Same as table 1, but the slope continues to y=2. 9 K ,0 The group velocities of the stable modes below the Rossby barrier are generally 2.5 times higher than for the FB model (see fig. 3.18), while the separation of the intercepts are magnified by 1.3. The group velocity curve for the fastest stable mode of the first kind has a flater profile over y2, offset by an enhanced curvature at small slopes. The group velocity curve for the second stable mode has lost its flat character for slopes bigger than unity; it wanes slightly. 48 None the less, its shape does not change qualitatively in a significant way. The separation between intercepts is marked by much higher curvature for small slopes. Fitting curves with simple functions with three parameters gave satisfactory results for the first type fastest mode: (>2+.1608) .2826arctan-1.475 and the intercept separation: .5991 arctan (ra+.1137) 1.021 However, the best fit for the second type mode was (/2-.03958) .2211tanh-which leaves significant residuals. Shelf-break with slope ends: a t f r o n t e n d a t y = 2 .5761 0.4 „ 0.2 O CL o oo 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 S e p a r a t i o n b e t w e e n s t a b l e m o d e s a t K = 0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Figure 3.18. Upper: Interpolated group velocities for a very long wave. Lower: Interpolated intercepts for the second type of stable mode. The absciss is the shelfbreak slope y2. Note that the scales for the two curves are different in the upper panel. Fig. 3.19. Sketch of the topography for a ridge of symmetric slopes and feet at y=-l and y=2, respectively. The changes for this model, with respect to the FB model, are that the slope continues beyond the front on the right-hand side and a slope occurs on the left-hand side for -1 < y < 0, with slope parameter y,. The choice of a ridge with equal slopes, say y = yx = y2, simplifies the discussion by reducing the parameter space. However, the symmetry does not influence particu-larly the dynamics of the system since the parameters (Z a) in the confluent hypergeometric equation are inferred separately for the two sides; the patterns are not symmetric in their real components. The most prominent changes are the apparition of stable modes above the upper Rossby barrier. For very long waves, the real component falls quickly and likely reaches zero when £=0. Inside the Rossby range, the sigmoidal stable branches are back combined with a symmetric pattern of marginal and parabolic branches and other stable curves with rightward direction in both marginal regions. At y=.25 (see fig. 3.20), a prominent change, compared to previous models, is a second kind of stable branch, which now has negative group velocity near £=0 and an intercept roughly between -.10 and zero. The positive Doppler-shifted nondimensional phase frequency reaches 0.064 at about AT=0.7. The unstable pattern is similar to that of the shelfbreak model with the abyssal plain at y = 2, with the exception that branch D is not rooted in the lower Rossby barrier. The end of branch D extends up to £=9.5, which allows to observe this branch carries the largest growth rates occuring similarly on branch A. T = W / R O +Ro -Ro 50 0.0* 0.03 0.02 0.01 |-•o.oo -0.01 -0.02 -O.03 -0.04 -0.05 Phase frequency +Ro 3 3  3 i 3 A 3 3 I 3 3 ' 5 5 5 *• 3 3 3 <? i 1 = 5 i 5 * 5 • • 5 3 ? S } 5 ° 4 4 E 2 3 , „ , „ „ _ , 2 2 2 2 -Ro 5 S 5 S 5 5 5 5 5 8 0.03 0.01 h Growth ra te Ro 3 ,A 3 3 B S 7 8 0 J ' , ' I 5 B 3jl5' 4 * 4 * • I ! ' 1 * E Figure 3.20.Same as fig. 3.5, but for a symmetric slope ridge with y=0.25 and feet at y=-l and y=2. At y=0.4 (see fig. 3.21), the stable branches frequencies outside the Rossby range increase by about 50%. The pattern of unstable branches is similar to the FB and shelfbreak models with the exceptions that: the tail of branch B is disturbed at A=2, causing a surge in its growth rate and a new locus, B', similar to the tail of branch E, but of opposite frequency; a peak in the growth rate occurs also in the tail of branch C, which suggests it was distorted and broken. 51 +Ro 0.04 0.03 -0.02 -0.0, --Ro -0.01 \--0.02 -0.03 h Phase f requency + Ro i > B • ' , 4 4 4 4 4 4 4 4 4 4 • o.oo k" 5 3 3 3 3 3 3 3 3 ., 3 3 3 iSu c „ 1 ! , n i ! ! ! ! E B . . 5 5 » 5 » « 5 S 5 c B 1 1 1 £ I I I 6 - R o T=( j /Ro -0.04 0 1 2 0 0 3 r Growth ro te Ro u 1 1 ' 5 1 4 6 1 8 5 8 III12 , 3 3 3 8' i 3 !32 9 K '° Figure 3.21. Same as fig. 3.5, but for a symmetric slopes ridge with 7=0.40 and feet at y=-l and y=2. At 7=0.63 (see fig. 3.22), two new stable branches appear above and below the Rossby barriers, respectively. Inside the Rossby range, a sigmoidal branch reaches the parabolic branches, but does not cause distortion. For the unstable pattern, the only differences with respect to the previous shelfbreak model are: the locus B', associated with branch E in the previous shelfbreak model, tends to coalesce with branch E at about A=2, rising their respective growth rate curves; the tail of branch E is extended smoothly to the origin; the previous shelf model showed a coalescence of branches E and B, but it appears as a crossing for this ridge model and causes the swap of labels between the tails with respect to the previous slope. The tail of branch C also has a growth rate peak at about £=4.5, characteristic of some critical point. For higher slopes, the number of stable modes increases on both regions outside the Rossby barriers and all stable branches move to higher wavenumbers and frequencies. The patterns of unstable branches for symmetric ridges are generally close to those for the previous shelfbreak models with equivalent slope. Meanwhile, their differences with respect to the FB models, are similar to those of the previous shelfbreak models. 52 Phase f requency +Ro j j ' ' ! i 33333 3 3333333J 1 A I 3 3 3 3 3 B 5  5 5 » 6 S E • 4 4 4 4 4 4 4 C 2 2 2 2 2 2 1 0 - Ro 0 1 2 3 0 0 3 f Growth rate Ro 1 1 * C , ' 332 1• S3,4' 2 " S 3 3 5 E A , 0 „ 3 * J * c 8 9 | ^ 10 Figure 3.22. M r _ — _ Same as fig. 3.5, but for a symmetric slope ridge with y=0.63 and feet at^=-l and>,=2. D 0 + f Phase frequency +R0 0.03 0.02 a o i -0.01 -0.02 h -0.03 —I L_ c o o ki; s«11 1 z I s s 5 s 5 B-8 8 6 I 3 3 3 3 3 4 4 4 C T = u / R o -0.0* 0 1 2 3 " • " f Growth ra te -Ro -0.751 / Ro A B Hi c ,2 t3' XL n 2 2 2 8 S 6 6 8 S S , 1 E ; I! 2 2 5 B" K -1 1 1 1 Figure 3.23. Same as fig. 3.5, but for a symmetric slope ridge withy=1.0 and feet at >»=-l and y=2. 53 At y=1.0 (see fig. 3.23), the only notable changes with respect to the shelfbreak model with foot at y=2, are a flater distribution of the growth rate peak at the tail of branch E and the reapparition, through one point, of branch D. Phase f requency +Ro 0.04 0.03 +Ro 0 0 2 0.01 I „™,L t 2 2 2 2 2 2 % 2 -0.01 -0.02 -0.03 , 5 E -0.0+ 0 0 3 ° Growth ro te x . - R o 3 4 Ro j A. B 2 2 2 2 2 3 , 2 3 ? ' B 5 3 2 2 2 5 . 5 3 3 5 i n i 2 2 2 , • 5 § E C 6 9 K 1 0 Figure 3.24. Same as fig. 3.5, but for a symmetric slope ridge with y=1.6 and feet at y=-l and y=2. At 7=1.6 (see fig. 3.24), the locus B' is significantly more unstable at its end (K=4), than for the previously labelled E' in the shelfbreak model with the same slope. The tail of branch E starts at a smaller wavenumber; its growth rate curve is similar to those of the B branch and B* locus and suggests it could originate from the origin despite a real part approching the lower Rossby barrier. The top growth rate on branch E is 0.011 at A=6. The growth rate curve for branch B peaks at AT=5.25 with a value of 0.01. Branches E and B get more stable at higher wavenumbers until they coalesce at £=8.5. The resulting branch has a sharply increasing growth rate, which converges with those of branches A and C. 54 0 0 *r Phase frequency 0.03 0.02 0.01 o.oo -0.01 -0.02 -0.03 i ! ! ! ? » . 2 1 11 • • . . . A; 3 , 2 2 2 1 2 5 5 5 5 s 5 B 3 , 4 4 E ' 'e 4 * 4 4 4 4 4 - R o i i i 4 9 6 I 1 7 8 9 Figure 3.25. Same as fig. 3.5, but for a symmetric slope ridge with y=2.5 and feet at,y=—1 andy=2. At y=2.5 (see fig. 3.25), the unstable pattern is again similar to the one from the shelfbreak model. The E branch is more detailed with a sharp peak at £=6 and is generally more stable than other branches for higher wavenumbers. Branch B still peaks at £=5.5, but undergoes another rise of its growth rate at about £=7.5 before it converges with branch E. A tiny locus, denoted E', is significantly more unstable than other branches at about £=4 with a growth rate of one third of the Rossby number. The locus A' extendsto the origin. 55 3.5. Shelfbreaks and symmetric slope ridges with long slopes 3.5.1 Shelfbreaks with foot at y=10 and Ri=57: With the foot of the slope in the farfield, much lengthier computations are needed for the left-hand boundary conditions and the workload increases with the wavenumber and as the frequency approaches either of the Rossby barriers. Therefore, the analyses will be limited to four values of y2: 0.25, 0.4, 1.0 and 2.5. The interior of the Rossby region provides a richer pattern of stable modes, which are better analyzed with the x-transformation. The points near the Rossby barriers are dense and difficult to calculate, because the series solutions are stiff and their differences are sensitive to any small change. Computations of the confluent hypergeometric solutions were quite lengthy and would have been improved by the direct use of a stepwise integration (such Runge-Kutta 4) with an adaptive stepsize control. Missed solutions are another side effect of the density of solution due to a lack of resolution on the initial grid. The points with |x|>0.7 are not plotted for clarity. The symmetric pattern of marginal and parabolic branches seen for shorter slope persists here, although it is evanescent for Y^IS. Other stable modes consist of non-overlapping sigmoids starting either from the origin or from the lower Rossby barrier to converge asymptotically on the upper Rossby barrier; and left-oriented parabolae for high wavenumbers. The left-oriented parabolae are not different from those observed for the shelfbreak model with shorter slopes but the sigmoidal stable branches, although confined in the same region, are densified by the longer slope. Outside the Rossby region, the primary effect of the longer slope on the stable modes is an enhancement of the group velocities at small wave numbers. The nonzero intercept of the last branch is moved to a much larger value. The stable branches along the lower Rossby barrier are stretched and undergo a reversal of their group velocity. Also, stable branches that could have been masked by the vicinity of the Rossby barrier are stretched to become apparent. For large wavenumbers, the stable modes remain nearly the same. The pattern of unstable branches is similar for this longer shelfbreak slope. 56 For /2=0.25 (see fig. 3.26), branch G observed for the shelfbreak at a shorter slope is absent. Branches B and C, which were found at shorter slope to have fluttering growth rates curves in the high wavenumbers, have points with two different growth rates for the same frequency; the fluttering part turns into two curves of growth rates and the highest ones have values near those of branches A and D. It can be noted that branch E, which is nearly imaginary, has a linear growth rate for wavenumbers where stable modes in the central area are sigmoids starting from the origin, a nearly constant growth rate for wavenumbers where the sigmoids originated from the lower Rossby barrier and a decrease in growth rate for wavenumbers in the range of transition from sigmoidal to parabolic stable branches 4.5<AT<5.5. The cutoff of branch E gives birth to the parabolic stable branch. a 0 4 r Phase frequency 0.00 r 0) -0.05 0.0 - M -0.4 -0.6 -OS Phase frequency of stable modes ;"; - _ (excluding inside Rossby region) 0 0 3 1 1 3 4 5 8 Stable frequency/ Rossby number -o .o i |--0.02 -0.03 h " -0 .04 -fRo / 0.00 K 1 1 ! s 1 2 l ' 1 1 > « > ' ' 1 1 I t E A 4 4 4 B „ ! I H 1 ! ' 1 ! c -Ro 0.03 Growth rate Ro 5 V 4 ^ § -K 9 K '° Figure 3.26. Shelfbreak with ifr=57, y2=0.25 and foot at y=10. Left upper panel: the stable modes outside the Rossby region. The dashed line is the lower Rossby barrier. Left lower panel: the x-transformation of the frequency for stable points inside the Rossby region. Right upper panel: the real part of the unstable points. Dashed lines are the upper and lower Rossby barriers. Right lower panel: the imaginary part of the unstable points. 57 For y2=0A (fig. 3.27), in contrast with the shorter slope, the decline of the growth rate for branch E can be observed. Its cutoff occurs at £=6.5. The unstable branches B and C have points with double growth rates for same frequency in the high wavenumbers. Inside the Rossby range, the transition from sigmoidal curve to a parabolic curve occurs at about £=5.25, which coincides with the transition of the unstable branch B into positive frequencies and with a peak in the growth rate of branch E. 0.0 to -0.1 -0.2 Phase frequency of stable modes (excluding inside Rossby region) 0.04 0.03 -0.02 -0.01 0.00 -0.01 -0.02 -0.03 Phase frequency +Ro i 3 3 3 3 3 3 • 3 3 3 3 ^ 4 4 * « 4 « « « « «B \ } ' V * \ \ 4 4 4 4 « ^ ' ' E C j 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 5 S 5 5 5 . - R o Growth rate / Ro 4 3 § - 5 • c 4 * 3 J « » 2 3 2 8 2 B , 3 I 1 * 1 1 ' 1% I 3 1 3 3 » 4 2 2 ? . E 1 K Figure 3.27. Same as fig. 3.26, but for a shelfbreak with y2~0A and foot at y=10. For /2=1.0 (fig. 3.28), in contrast with the shorter slope, the tail of branch E starts smoothly from the origin and no peak of instability occurs near £=3.5. One point signals a return of branch D for very high wavenumbers. The transition from sigmoidal curves to parabolic curves occurs near K=6. 0.0 -0.1 - 0 2 1 -OS 0.6 0.4 h 0.0 - O J h -0.4 -0.61-- 0 J 58 Phase frequency of stable modes (excluding inside Rossby region) 1 2 3 4 5 6 Stable frequency/ Rossby number 0.04 0.03 0.02 0.01 Phase frequency +Ro o.oo < ' 5 -o.oi --0.02 -0.03 5 5 5 5 5 5 5 5 5 5 g ' l i l t \ - R o 1—-0.04 K 0.03 1 2 Growth rate 3 4 5 / Ro A B 4 4 4 2 1 » 1 1 5 K 9 K 1 0 Figure 3.28. Same as fig. 3.26, but for a shelfbreak with y2-l.O and foot at y=10. For y2=2.S (fig. 3.29), in contrast with the shorter slope, branches A and B span smoothly from the origin through the wavenumber range. The growth rate of branch A is slightly depressed about £=9.5. Branch E is distorted twice with peaks of growth rate at £=5 and £-6.25. 59 0.0 r . 0) 0.6 • 0.4 • 0J • ajo • ^hase.frequency of stable modes (excluding inside Rossby region) - i — ' — r -2 3 Stable frequency/ Rossby number ° ' w r Phase frequency 0.03 h 0.02 0.01 0.00 -0.01 --0.02 --0.03 - i 1 1 ' 1 < r-5 6 7 8 -0.04 0.03 0.02 +R.0 li 3 J 3 3 > i 3 3  3 3 3 3 3 2 2 2 2 2 J 2 2 2 2 2 2 2 2 2 2 2 ' 1 1 1 1 . 1 1 A 3 3' 2 B 2 -Ro I —I I I 1 1 1 1 T 1 1 2 3 4 5 6 7 8 9 10 - i ' 1—1 1 1 — r -Growth ra te / Ro a 3 3 3 ' 2 2 2 1 3 A 3 B - i 1 1 1 Figure 3.29. Same as fig. 3.26, but for a shelfbreak with y2=2.5 and foot at y=10. 60 3.5.2 Symmetric slope ridge with R/=57, feet at y=-9 and y=10: Farfield boundary conditions occur on both sides of the front, which cause the longest computations for the dispersion function. The analysis will be limited to y=T (see figs. 3.30 and 3.31). Stable modes occur on both outer regions. As the longer slope increases their group velocities for small wave numbers, it increases their density near the Rossby barrier and deforms them in the same manner as for the shelfbreak case by bringing stable mode extrema to smaller wavenumbers. The top stable mode with positive frequency is observed with a finite intercept extrapolated to 0.56. Exploration for very small wavenumber reveals it actually peaks at 0.5232 with a flat region between 0.\<K<0.14, and a tail attached to the origin. Similarly, the top stable mode with negative frequency peaks at -0.5357 at £=0.11. For the stable points inside the Rossby range, the usual symmetric stable pattern is masked by a less orderly but dense skew pattern. The structure for high wavenumbers (£>7) and moderate frequencies is relatively unaffected by the length of the slopes. Apart from a minor kink on branch E at £=3.25, which causes a peak of 0.0095 in the growth rate, no major distortion occurs similar to the shorter slope. The unstable pattern is very similar to that of the shelfbreak with the longer slope: the dispartures being the distortion outlined previously for shorter slopes and the absence of one point of branch D. Figure 3.30. Symmetric ridge with Ri=57, y=1.0 and feet at y=10 and y=-9. Upper panel: the stable modes outside the Rossby region. Lower panel: the t-transformation of the frequency for stable points inside the Rossby region. 0.04 r 0.03 0.02 h 0.01 O.OO h -0.01 -0.02 -0.03 \-P h a s e f r e q u e n c y «. 4 4 * * * * • * 3 3 3 3 3 3 3 3 3 3 + R o A 4 • 4 * 4 * 3 3 3 3 3 3 3 3 3 3 ? 1 1 i 1 1 , i , , E 1 1 ' 1 c z z -0.04 O 0.03 r 0.02 h Q.01 1 I — 1 2 G r o w t h r a t e - R o j— s e Ro 1 3 3 2 B A 4 4c 4- 3 4 3 1 8 3 1 t * 1 1 1 3 4 ' ~ T — 8 7 B S 10 K — i 10 62 Figure 3.31. Symmetric ridge y=1.0 and feet at y=10 and y=-9. Right upper panel: the real part of the unstable points. Dashed lines are the upper and lower Rossby barriers. Right lower panel: the imaginary part of the unstable points. CHAPTER 4 63 THE LOMONOSOV RIDGE 4.1 The observed front during ODEN 91 expedition: Reality posed several difficulties for the application of the model. From the observations of temperature and salinity during the ODEN 91 expedition (Anderson et al. [1994]), the various stratifications, the surface salinity gradient and the presence of a warm poleward subsurface current over the ridge crest cause several bows to the interface between the Eurasian and Cana-dian water masses. The ridge itself presents a blunt crest that can be represented as a 35 km wide plateau. It was observed from the isotherms from the ODEN 91 expedition (fig. 4.1) that the warm core of the poleward current, at a depth of 300 decibars with a temperature maximum higher than 1°C, is above the upper part of the ridge slope on the side of the Amundsen Basin. An abrupt diving of the isotherms signals a front attached to the upper part of the slope on the Canadian side of the Lomonosov Ridge. The strong meridional temperature gradient between 0 and 300 km is the Atlantic warm layer flowing along the southern boundary (Eurasian Shelf) over the Nansen Basin. Figure 4.2 gives details of the top 1000 decibars. The Gakkel Ridge induces a front between 100 and 300 km that prevents poleward circulation of the Atlantic water below 300 decibars, but allows some lateral intrusions between 200 and 300 decibars. It can be interpreted from the contours that, as the warm Atlantic layer penetrates deeper along the southern boundary (across the transect), it crosses the Gakkel Ridge to pass into the Amundsen Basin. There, the intrusions from the Atlantic layer recirculate freely toward the pole. The intrusive layer at the depth of 300 decibars benefits the most from lateral recirculation, and gives the temperature maximum observed above the Lomonosov Ridge with an elongated tail over the Amundsen Basin. Remnants of a second major fossile intrusive layer that occurred at a depth of 220 decibars can be perceived above the Lomonosov Ridge as a deformation of the 0.7°C, 0.8°C and 0.9°C-contours above the temperature maximum. 64 Potential temperature contours (°C) 1—' • — - i r • 1 1 1—• r 1 r -0 100 200 300 400 500 600 700 Southern boundary Distance (km) (g) Poleward Figure 4.1 Contour plots of potential temperature observed during the ODEN 91. The dashed line from the top of the Lomonosov Ridge shows the location of the front beneath the pycnocline. The density profile computed from potential temperature and salinity observations at station 23 of ODEN 91 (personal communication from Peter Jones, Bedford Institute of Oceanography, Canada), reveals that the water column between 236 and 261 decibars was unstable and in the process of vertical mixing (see fig 4.3). The vanishing of the second major fossile intrusion can be explained by attacks of new intrusive layers at the microscale. A consequence is that sinking water is replaced laterally by fresher and colder water of Canadian origin that forces its way through the front between 265 and 300 decibars. The process can be perceived in the subsurface temperature contour plot in which the 0.5°C-isotherms bordering the warm core do not quite meet near 300 decibars, but penetrate through the front and plunge below 450 decibars, 100 km away into the Canadian basin. 65 Potential temperature contours (° C) Pressure Station number 20 24 (decibars) 9 10 11 12 13 14 15 16 17 18 19 21 22 23 25 26 Figure 4.2 Contour plot of potential temperature observed during the ODEN 91 in the top 1000 m. The near vertical dashed lines show vortex sheets in the neighborhood of the Lomonosov Ridge. The vortex sheet on the right is the frontal interface beneath the pycnocline. The arrows show a sinking of water of Atlantic origin into the Canadian water assembly. Pressure (dbore) T (dog. Celsius) 0 . 6 0 0 . 6 5 0 . 7 O 0 . 7 5 0 . 8 0 O . B 5 0 . 9 0 .1 1 3 4 . 7 8 3 4 . 8 0 3 4 . 8 2 3 4 . 8 4 3 4 . 8 6 3 4 . 8 8 SljJ.S.U.) Figure 4.3 Profile of potential temperature (T), salinity (S) and density-1000 (o) between depths of220 decibars and 320 decibars from station 23 of ODEN 91. Arrows indicate levels where vertical gradient of density is positive, which triggers vertical mixing. The water column between 236 and 262 decibars is going to undergo vertical mixing. 66 Beneath the warm core, deformations of the isotherms suggest a vertical vortex sheet, standing along the Eurasian edge of the Lomonosov Ridge, that is distinct from the vortex sheet associated with the front on the Canadian edge. Albeit it constitutes a front, it stays a minor feature since it does not separate two distinct water masses. The poleward current is associated with positive temperature anomalies (downward shift of the contours below 300 dbars) above the slope, however the combined drag due to friction with the blunt crest and crossfrontal momentum diffusion (from microscale turbulence) generates negative temperature anomalies at the boundary of the Eurasian water mass, directly above the ridge. In the case of a narrow crest, both vortex sheets would join to enhance the momentum difference across the front, which is a destabilising factor for a front, and the warm core would be closer to the main front and exposed to greater erosion from microscale intrusions. In the observed case, the stability of the front benefits from the buffer range made by the drag of the blunt crest. However the addition of a secondary front is not justified in the model because the ground friction and turbulent viscosity associated with very high wavenumbers are not taken into account. Another striking feature is that water temperature above the thermocline varies differently than in the layer immediately below. While the surface conditions in an open sea depend mainly on atmospheric features, the icepack in the Arctic imposes a boundary condition: the sea surface temperature must be at the freezing point. The freezing point varies weakly with the surface salinity. Surface water is fresher on the Canadian side: it is lighter than Eurasian surface water. The outcropping front bends toward the Eurasian side, causing a poleward boundary flow of Canadian water. Albeit the vertical gradient of temperature is generally small in the surface mixed layer under the icepack, the poleward boundary flow causes sufficient warm advection to bring the -1.7°C-isothermal to the surface. 67 4.2 The modeled ridge with a blunt crest: Among the major simplifications of the model, layer stratification is not included and the model does not take into account horizontal density variations inside the layers. The front slope depends on the horizontal pressure gradient across the front. Tremendous horizontal density variations occur mainly above the pycnocline, but the difference between vertical gradients of density contributes at greater depths. If the stratification difference could offset horizontal density variations (i.e. if the total mass above does not vary), the front would become vertical at that level. That level could be qualified as barotropic. There are as many barotropic levels as there are turning points in the frontal interface. Beneath the pycnocline, vertical gradients of density become nearly constant. If internal waves over the pycnocline are not considered (rigid lid approximation), the proposed model can be applied by replacing the sea surface with the deepest barotropic level. The curvature of the frontal interface will be neglected. The blunt crest is modeled as a flat top. The modified pressure perturbation field and its derivative field, the left boundary conditions for the frontal system, are: Canadian Edge 0" ~eKA e-KA - "1 0 " (PA K_ eKA -e~K& 0 Eurasian Edge where A is the width of the plateau in units of the frontal length. The lowest barotropic level was estimated to be about 600 m beneath the surface. H is redefined as the distance between the crest and the barotropic level, that is roughly 750 m. From the contour plot, the estimated frontal length was 15 km. Thus the geometric ratio, L/H, is 0.02 km/m. The plateau is about 35 km wide or 2.333 units of frontal length. The floors of the neighboring basins are considered flat and extended to infinity. The floor of the Amundsen Basin is estimated to be 3100 m lower than the plateau. The foot of the ridge on the Eurasian side is about 75 km away from the edge of the plateau (5 units of frontal length). The slope times the geometric ratio gives yx =0.827. On the Canadian side, the floor of the Makarov Basin is 2800 m lower than the plateau. The estimated foot of the ridge is 50 km (3.333 frontal length 68 units) away from the edge of the plateau, which gives y2=l.\2. The abyssal plain starts 2.333 units of frontal length away from the end of the front. I. f=0.000146 rad/sec Eurasian side Slope=41.33 km/m y 1 = 0.827 Amundsen B a s i H=750 m Abyssal p l a i n 600 m Ba r o t r o p i c l e v e l *i_=0.02 kiri/m H Canadian side Slope=56 kifi/m T2 =1.12 Makarov Basin Abyssal p l a i n 75 km 35 km 115 km_ 35 km y=-7.333 y=-2.333 y=0 y=l y=3.333 Figure 4.4 The model representation of the Lomonosov Ridge frontal system. The flow speed being unknown, the Richardson number was estimated from Ri = (fL)2 / g' H = {p/Ap)x(jZ)21gH = Q.671kgm~31Ap, with p=1028.19 kg/m3 as the density scale below 800 decibars. Ap should be an average density difference below the barotropic level between the Eurasian and Canadian water assemblies. It was estimated from temperature and salinity observations at 1000 decibars at stations 23 and 26 that A/?=0.0102 kg/m3 and #/=66. The search of eigenvalues in the complex plane leads to only one unstable mode (see fig. 4.5), varying linearly with the Rossby number for small and moderate wave numbers: = Ro f-Q. 5466^ 1 'i J + 0(Ro2) for K<3 or #o<0.0227. 0.2907 The group velocity of unstable long waves is |vjJ|=.004141/Z,= 7835 km/day. Behaviour changes after £=5, with a Doppler-shifted frequency oscillating between -.01435 and -.015585, then stops before K=7. The growth rate also oscillates between .0096921 and .008145, but it does not fit a continuous curve beyond A=5. The inability to find eigenvalues, at the limit of their accuracy in single precision, that satisfy the conditions with high accuracy (i.e. low score) 69 at high wavenumbers may explain the spread. However this is not always the case, the score is about 10"4 for A=6, and the eigenvalue is (-.0146538, .009692), which is a dominant eigenvalue. That means a frontal wave with a length of 2387 km would grow at a rate of 13% per day. The least growth rate for A>5 is 10.82% per day. The group velocity is positive in the range 5.24<£<5.9, with a value of 0.2323 km/day, and beyond K=6.5 with a value of 0.4181 km/day. Unstable mode for Lomonosov Ridge with plateau 0 1 2 3 4 5 8 7 8 wave number K Figure 4.5 The real and imaginary parts of the eigenvalue for the unstable mode in the Lomonosov Ridge case with flattened crest. Stable modes are present both inside and outside of the Rossby range. The pattern of stable modes outside the Rossby range (see fig. 4.6) is, as expected, most similar to the one for the case of a symmetric ridge with feet placed nine frontal lengths away from the front and slope equal to unity (fig. 3.30). As it could be expected from feet placed closer to the front ends, fewer branches are observed and their density near the Rossby barriers is decreased. Also fewer branches have points of zero velocity group. However, the higher Ri means the branches can expand to larger frequencies uniformly over the wavenumber range. There is a much higher density of modes beneath the lower Rossby barrier than above the upper Rossby barrier. The prominent and unexpected feature is the top positive branch with an apparent finite intercept at 0.6333. Explorations for small wavenumbers in the order of 10"2 did not resolve any peak. The top negative branch peaks at -0.5674 at about K=0.18. 70 0.6 » 0.4 cr e £ XL tn 0.2 k. J D & 0.1 |-o O 0.0 Stable modes f o r the Lomonosov Ridge (excluding modes inside the Rossby range) Ro -Ro 6 7 8 wave number K 10 Figure 4.6 The stable modes outside the Rossby range for the Lomonosov Ridge with a flattened crest. The interior of the Rossby range provides a richer pattern of stable points (fig. 4.7). The regions 0.75 RO<\G>\<RO were removed from the figure because of incompleteness and lack of clarity. Part of the stable pattern is still reminiscent of the symmetric pattern in spite of changes in the slopes length (marginal and parabolic stable branches), even if the slopes do not have here the same value and lengths. The parabolic branch starts at &=5, which coincides with the onset of the oscillating pattern for the unstable mode. Among the stable branches that depend on the ridge feet locations are two families of sigmoidal curves which converge asymptotically toward the upper and lower Rossby barrier, respectively, after starting either from K=0 or from the opposite Rossby barrier. Each family overlaps the symmetric pattern or the other family, but no curve overlaps another curve from the same family. Among the sigmoidal curves converging toward the upper Rossby barrier, a transition occurs beyond K=7 to parabolically shaped curves. 0.75 0.50 0.2S 0.00 r 71 T = = W / R O -0.25 -o.5o \=- =-= - - =-.;-;: - - --0.75 h 0 1 2 3 4 5 6 7 8 9 10 wave n u m b e r K Figure 4.7 The stable eigenvalues with the transformation inside the Rossby range for the Lomonosov Ridge with a flattened crest. The regions 0.75<|x|<l. were avoided for clarity. 72 4.3 The effect of a sharper crest on the Lomonosov Ridge: The model is set back to a triangular ridge by setting the plateau width to zero. For the stable modes outside of the Rossby range (see fig. 4.8), the change effects only the small wave numbers, mainly by suppressing the finite intercept of the fastest positive mode, which therefore peaks at 0.28 at about £=0.45. However, it effects greatly the modes inside the Rossby range (see fig. 4.9) with the substitution of sigmoidal curves crossing the region |x|<0.3 and K<6 by C -shaped curves resulting from the merging of the two types of sigmoids. The regions of high wavenumber and upper marginal branches are relatively uneffected. Inside the range 2.5<£<5, is a zone of broken C-shaped curves: the group velocities are too large to allow continuous curves. Stable modes for a Lomonosov Ridqe with narrow crest n . . . . . . . (outside Rossby range) 0.2 ' ' 0.1 0.0 -0.1 -0.2 -0.3 -0.4 -0.5 Ro 1 ' l| . ' ' 1 I I M I I I I ! I 1 1 I I I I I I I 1 -Ro , i i i 0 1 2 3 4 5 6 Figure 4.8. Stable modes outside the Rossby range for a narrow Lomonosov Ridge. 73 r-co/Ro 0.75 0.50 0.25 0.00 -0.25 -0.50 -0.75 h-0 1 2 3 4 5 6 7 8 9 wave number K Figure 4.9. Stable modes with the t-transformation inside the Rossby range for a sharp Lomonosov Ridge. The most important change, from the plateau removal, is the splitting of the unstable mode into several branches (A, B, C, D and E on fig. 4.10) with a pattern similar to the previous results for shelfbreaks and symmetric ridges. Branch A is the most unstable beyond A=6 with growth rates about 0.0105 or a growth factor of 14.2% per day and another likely growth rate curve about 0.0125 or a growth factor of 17% per day for K>7.5. The frequency of branch A is always near Ro/4. Both branches A and B have unsmooth growth rates that are double for some wavenumbers. Branch C is much smoother than A and B, and the part for K<3 has frequency equal to -Roll, which is in the region of the lower marginal stable branches. The point D is also placed in the region of lower marginal stable branches. This result is in opposition to the previous observation that the frequencies of unstable branches are confined outside the region spanned by the marginal stable branches. Branch E is purely imaginary and has the lowest growth rate. 74 0.02 a) 0.01 -0 .02 N a r r o w c r e s t m o d e l of L o m o n o s o v Ridge Dopp le r s h i f t e d f r e q u e n c y / R o s s b y n u m b e r -3 Sf 33 33 ^  3 3 3 3 3 3 "£"333 3 3 . 3 * 3 3 - 3 3 " W 4 :'V" 2 2 1 1 1 1 1 1 1 1 1 1 1 1 E * ~4 ... . * 2 2 2 2 2 2 2 2 2 2 2 22 2 2 2 22 2 2 2 2 2 2 B < 4 4 4 4 4 4 4 * 4 * 4 C - R o / 2 b) 0.012 0 1 ,2 Growth ra te 0 .004 0.000 7 8 9 wave n u m b e r K R o / 3 3 j 3 3 3 a s 3 a 3 * . 2 * , 2 4 / 3 * 3 » 2 4 I 2 , 4 a 2 1 « 3 f , 1 4 2 1 1 1 1 1 1 1 E 10 Figure 4.10. The unstable modes found for the sharp crest model o f the Lomonosov Ridge: a) the Doppler-shifted frequency divided by the Rossby number, b) the growth rate. CHAPTER 5 DISCUSSION AND CONCLUSION 75 Using power series eigensolutions instead of a stepwise integration of finite difference equations, the model faithfully reproduces Orlanski's [1968] results (flat topography case with i?/=3) for small and moderate wave numbers. However, it does not find instabilities beyond what we define here as K=A or #o=0.7. The discrepancy can be explained by Orlanski [1968] who mentioned a poor convergence of his technique in that region. The method used here is quite limited beyond AMO. At high Richardson number (/?/=57) for a flat topography configuration, the model shows symmetric unstable branches and a symmetric pattern of stable branches within the Rossby range. The most unstable branches have frequencies midway between the axis Re((o)==0 and the Rossby barriers. The same branches bear the highest growth rate at any wave number; their growth rate curve increases steadily in the observed range, but extrapolation from the curvature suggests a maximum growth rate of about 0.03 beyond the observed range. The introduction of a slope beneath the front (FB model) creates stable branches of backward-propagating modes with nondimensional Doppler-shifted frequency greater than the Rossby number. The fastest among them has a forward group velocity at small wavenumbers with a finite nonzero frequency intercept. An increase in the slope moves the stable modes away from the lower Rossby barrier over all the observed range of wavenumbers. Another effect is the occurence of new stable branches along the lower Rossby barrier. These branches might have been there for lower slopes, but the roots would have been either too near the pole or rejected because the solutions overflowed or were not accurate enough in the differential equations or their separation was not enough for observation by the bisection technique. The greater slope unmasks them by increasing their distance from the pole. An increase in the length of the slope has the same effect for small wavenumbers. The smaller the wave number, the greater is the effect. Therefore, group velocities for small wavenumber are enhanced by a longer slope. For a long slope, the stable modes get denser along the lower Rossby barrier. It can be extrapolated 76 that the number of stable modes for an unbounded shelfbreak slope becomes itself unbounded and the majority of the roots would concentrate in a domain of poor convergence. The use of asymptotic solutions for an unbounded slope with Tricomi function was rejected because computations would be too tedious for the wanted accuracy and the program was fated to crash when the storage capacity of arrays would be overwhelmed by the cumulated roots. This problem was not fixed since Nature itself does not support such a case. Forward-propagating roots were found during the exploration of the FB models at the very limit of \+Ro, for moderate wavenumbers and large slopes. The locus lies along the line on which the model becomes invalid. For larger slopes, the locus spreads in its extent and undergoes duplication of its points, followed by a densification at higher slopes. Those points at the limit of validity have been dismissed for convenience. However, they deserve further investigation for the steep and long slopes combination. No other stable mode with subgeostrophic positive frequency greater than Ro occurs for the FB model or the shelfbreak model. This kind of stable mode occurs for a ridge model, although the relative accuracy of roots is much poorer than for the backward-propagating modes, in particular for high wavenumbers. The length of the slope influences the group velocities for the large wavelengths. For the stable modes at high Richardson number, the persisting symmetric pattern inside the Rossby range consists of branches running nearly parallel to the Rossby barriers, which start from the origin and have values above Ro/4 (marginal branches), and of a pair of branches starting at a high wavenumber and forming a parabolic shape (under the t-transformation). The marginal branches appeared in a countable manner for low Richardson number. For greater Richardson number, they become more numerous and denser near the Rossby barriers. In the neighborhood of the poles, the series solutions behave closely and form a stiff system where numeric truncation errors remove significance of close roots, which makes the set of marginal branches uncountable for high Richardson number. Although the parabolic branch is observed for the flat model at high Richardson number, it is not found for low Richardson number. 77 Instead, a branch of mixed instability type, associated with Orlanski's B surface, is parabolically shaped and terminates with the birth of uncentred, parabolic stable branches. The symmetric pattern depends upon the slope beneath the front and the Richardson number. Its existence is not modified by topographical change beyond the slope. As the series solutions are built, later terms in the series become dominant for frequencies approaching the Rossby number or high wavenumbers (which slow down the convergence), and the memory of the boundary conditions (ratio of pressure field and its derivative at the front end) is deleted by numeric truncation errors. The introduction of a plateau for the Lomonosov Ridge did not change the symmetric stable pattern. However, the boundary conditions create asymetric stable branches that interact with the symmetric pattern. Distortions occur at the crossings between the branches of the persisting symmetric pattern and those of the pattern depending upon the boundary conditions: branches avoid crossing by redirection and path exchanges. This sort of scattering is accentuated with strong slopes. The unstable pattern remains contained inside the semicircle of radius Ro, in the complex plane. Its' features depend strongly on the slope beneath the front, but the unstable pattern is nearly same for the FB model, the shelfbreak with the foot in the near-field and the symmetric slope ridge with feet in the near-field. The largest observed growth rates were under 0.025 at about the end of the observed range. That is 17% below the largest values observed for the flat bottom case. Displacing the foot can create distortions for moderate wavenumbers (breaking a branch tail or causing its duplication), accompanied with a peak in growth rate. For longer slopes, fluttering growth rates occur and some unstable frequencies have two distinct growth rates. The introduction of a plateau between the slope and the beginning of the front in the Lomonosov Ridge model is a major change in the boundary conditions that reduces the unstable pattern into an unique branch. Gawarkiewicz [1991] found for his shelfbreak model, three main unstable features for Ro between 0.5 and 1.8 (57<£<200), Ri=57 and y2=6, with a top growth rate of 0.015. His main growth rate curves are probably the extension of the unstable branches A, B and C in this dissertation. To get these curves for very high wavenumbers and unbounded 78 slopes, he used in conjunction with the Hoskin's [1975] quasigeostrophic momentum approximation for the perturbation, an asymptotic expansion for the Tricomi function in the region of the unbounded slope. The unstable features were found by this simpler approach, but the accuracy can be questioned. However, he found that a change from the false bottom conditions to an unbounded slope condition causes little change in the most unstable features, which agrees with the observations for high wavenumbers when the foot is displaced. He also found that the most unstable feature has the most negative phase velocity. Here, the results are more mitigated: Branch A bears in several cases the most unstable eigenvalues with positive frequencies. Gawarkiewicz's [1991] curves are not as consolidated as those in this work. The application of the model to the Arctic marine front near to the Lomonosov Ridge is hindered by the lack of knowledge of the frontal length, L. Stratification is a major factor that was overlooked. However, the model was applied for that part of the water assembly beneath the pycnocline, for which the frontal length could be estimated and vertical variations are reduced, and assuming the surface variations are cancelled at the level the front becomes barotropic. The main finding is the reduction of the unstable pattern to one branch that terminates abruptly after K=7 (wavelength shorter than 2000 km). It provides growth factor of at most 13% per day for wavelength about 2400 km. BIBLIOGRAPHY 79 Anderson, L.G., G. Bjork, O. Holby, E.P. Jones, G. Kattner, K.P.Koltermann, B. Liljeblad, R. Lindegren, B.Rudels and J. Swift, Water masses and circulation in the Eurasian Basin: Results from the Oden 91 expedition. Journal of Geophysical Research, 99:3273-3283,1994. Barth, J. A., Stability of a coastal upwelling front. Part 1: Model development and a stability theorem. Journal of Geophysical Research, 94:10844-10856, 1989. Bjerknes, J., On the structure of moving cyclones. Geofys. Publikasjioner, 1(2), 1919. Boyce, W.E., and R.C. DiPrima, Elementary differential equations and boundary value problems. Wiley, fifth edition, 1992. Charney, J., The dynamics of long waves in a baroclinic westerly current. J. Meteor., 4:135-162,1947. Duxbury, A.C., An investigation of stable waves along a velocity shear boundary in a two-layer sea with a geostrophic flow regime. Journal of Marine Research, 21:246-283,1963. Erdelyi, A. (editor), Tables of Integral Transforms, volume 1. McGraw-Hill, 1954. Fj<|>rtoft, R., Some results concerning the distribution and total amount of kinetic energy in the atmosphere. In The Atmosphere and Sea in Motion. Rockfeller Institute and Oxford University, New York, 1951. 509pp. Flagg, C. and R.C. Beardsley, On the stability of the shelf-break/slope-water front south of New England. Journal of Geophysical Research, 83:4623-4631, 1978. Gawarkiewicz, G., Linear stability models of shelfbreak fronts. Journal of Physical Oceanography, 21(4):471-488,1991. Hoskins, B.J., The geostrophic momentum approximation and semigeostrophic equations. Journal of Atmospheric Sciences, 32:233-242, 1975. Kotschin, N. , Uber die stabilitat von margulesschen diskontinuitatflachen. Beitrdge zur Physik derfret Atmosphare, 18:129-164, 1932. McLaughlin, F., E. Carmack, R.W. McDonald, J.K.B. Bishop, Physical and geochemical properties across the Atlantic/Pacific water mass front in the southern Canadian Basin. Journal of Geophysical Research, 101:1183-1197,1996. Orlanski, I., Instability of frontal waves. Journal of Atmospheric Sciences, 25:178-200,1968. 80 Pedlosky, J., The stability of currents in the atmosphere and the ocean. Journal of Atmospheric Sciences,2\:20\-2\9,1964. Phillips, N. , Energy transformations and meridional circulations associated with sample baroclinic waves in a two-layer model. Tellus, 6:273-280,1950. Rainville, E., Intermediate differential equations. MacMillan, second edition, 1964. Rudels, B., E.P. Jones, L.G. Anderson and G. Kattner, On the intermediate depth waters of the Arctic Ocean, in The Polar Oceans and their Role in Shaping the Global Environment: The Nansen Centennial Volume, Geophys. Monogr. Ser., 85, edited by O.M. Johannessen, R.D. Muench and J.E. Overland, AGU, Washington, D.C. 1994. Solberg, H., Integrationen der atmospharischen stflrungsgleichungen. Geofys. Publikasjioner, 5(9), 1928. 81 APPENDIX A THE TRICOMI SOLUTION OF THE CONFLUENT HYPERGEOMETRIC EQUATION In this section, the solution of the confluent hypergeometric equation (1.17) will be discussed. Let us simplify the notation by using (z, a) instead of (zx,ax): (A.1) zpzz+(l-z)pz-ap = 0. The classical solutions can be constructed by series expansions about the singular point at z - 0 (or y=\JY\\ with the Frobenius method (see Boyce and diPrima [1992], section 5.6). The first characteristic solution assumes a regular singular condition at z - 0. Let n=0 Substitution into the confluent hypergeometric equation leads to r(r - l)a0zr~l + ra0zr~l + 0(zr). The balance of the coefficients at the lowest power of z leads to the indicial equation r2 =0. The multiplicity of the root indicates that only one solution has regular behaviour at the singular point, while the second kind has a logarithmic singularity. Since the singular point is not contained within the region 7, the irregular singular behaviour is physically possible. The indicial equation shifted by the power n is an+i(" +!)" + (" + ! ) - na„ -aan=0. Thus, the recurrence relation is n a„^ = ~an = J ~ | - — - r t f r . = Let a0 = 1. The first kind of solution is which is a confluent hypergeometric solution known as the Kummer function (see Erdelyi 1954]). The Kummer function can grow in an exponential-like manner for large z. n+a _ -pr i+a a0 F(a+n + l) 4"+l = (n + lf a" = ! i ( 7 ^ r ° = (n + lf T(a) 82 The ratio test: Urn n + l a z n 2 Urn n + a = 0, \(n + \)2 thus the radius of convergence of the solution is infinity. A special case of note is when the parameter a is a nonpositive integer, say -n. The ordinary differential equation becomes zpzz+(l-z)pz + np=0. which is the Laguerre's equation and the Kummer function becomes a finite series identical to the Laguerre polynomial Ln(z) of degree n: xFx(-n;l;z) m L„(z) = 1 + £ f[( 1 - £ ± I U (n > 0) = 1 (n = 0). The second kind of solution, which expresses a logarithmic singularity at z=0, is allowed in the system because the singular point is not part of region /. As the indicial roots are zero, the Frobenius method asserts that Pn(z) = Pi(z)\ogz + g(z), where Then PnAz) = PiAz){o^z+Pi(zVz+sz(zl Pu,zz(z) = Pi^(z)^z + 2pIyZ(z)/z-Pl(z)/z2 +gzz(z). The differential equation becomes 2Pi,z ~Pi/z + 8zz +(\-z)Pilz + Q-z)gz - &g = 0 or L. H. S,=ga + (1 - z)g2 - ag = pj - 2pj z =1F1(a;l;z)-2a lFl(a+ l ;2 ;z) = R. H.S. The power expansions of each side are: 83 r(a + n) z" o r ( a ) ( « ! ) 2 •21 n=0 r(g+n) z"-1 F(a) n ! ( w - l ) ! ^«!r ( a ){ «! (w + 1)! J ^ ( » ! ) 2 r ( a ) l n + l j - l - 2 « + £ ; r < » + 1 * " =! («0 2 r(a) for a not a nonpositive integer, or 1-2 H + l for a a nonpositive integer; Z,//.S.= zg" + g ' - (zg ' + « g ) = ft + X fe»+t + + g„(" + a)}zn. «=i A term-by-term comparison gives g„+i(« +1) 2 -&,(" + «) = 1-2 n + a n + l n i + a . . 5- for n > 1. That is _ n + a-l Sn ~ 8n-l 2 ^ 1-2 1 + a-l n 1) Therefore, the general solution of (A. 1) is (A.2) p(z)=lFl(aXz)(C\ogz + D) + cf^gn(a)zn, where C and D are constants of integration. Its derivative is pt(z)=alFl(a+l;2;z)(C\ogz + D) + C ^ ( « ; l ; z ) / z + £ « g „ z " + 1 =o(Clogz + D) a~\\z (A.3) c + — z i + IJi n = l V i = l ^ I J I J 84 For a particular solution that decays in the far-field, a ratio D/C must be found in order that p(z) is proportional to the Tricomi function (Tricomi [1927], Erdelyi [1954]): where XFX(a, 1;z)log z + £ ^ ( " + w ) * { y<W + a) - 2 y<» +1)} Ho («!) r (a) yKx)=-r \ogY(x). ax It suffices to set D/C = yA^cx) ~ 2 if/(l), that is ti n(n-l) where y ^ O . 5772156649 is the Euler-Mascheroni constant. Unfortunately, the Frobenius' fundamental solutions can grow rapidly and have similar asymptotic behaviour for moderately large z, and the numerical determination of the general solution becomes inaccurate. In the case where the slope stretches to infinity, the asymptotic behaviour of the solution is better observed in the far-field (large z domain) with the transformation rf=l/z. The ordinary differential equation becomes (A.4) rfpm + {T]2+ ij)Pri -ap = 0. Assume that the function remains smooth and does not have a logarithmic singularity at infinity. Then n=0 Substituting into the ordinary differential equation gives £an (n + s)(n + s -1) rf^ an (n + s) rf^X + 1 a„ (n + s) tf+* - a £ a n rf+s = 0. n=0 n=0 n=0 «=0 That is a0 (s -a)Tjs + rjsJ^ [an_x (n + s-1)2 +a„(n + s- a)]rf = 0. To be independent of TJ, s=a and an = - a , ^ = ^ = a o f l ( - 0 ( l + ( « - l ) / / ) 2 -if Thus, the far-field solution is 85 (A.5) P = -Jt \ + YXl(\ + {a-\)li)\-ilz) 71=1 1=1 The Tricomi solution decays in the far-field only when Re(a)>0. Otherwise, the asymptotic behaviour of the confluent hypergeometric solution has polynomial growth in the far-field. However, the overall solution is exponentially-damped, which fulfills the physical boundary conditions, whatever the value of a. In the specific case where a is a nonpositive integer, the asymptotic solution is proportional to the Laguerre polynomial of degree -a . The derivative of the asymptotic solution is A - ^ t + 0 ( , * . ) . For a more detailed expression, let q=pz and the differentiated ordinary differential equation is zqa + (2 - z)qz - (a + \)q = 0. With the transformation JJ = it becomes (A.6) ifqm+rjqv-{a+\)q = Q. Substitute q - Xn=0 ^» m t n e ordinary differential equation: 0 = ^b^n+r)(n+r-l)7f+r+l+J^bSn+r-a-l)rf+r n=0 «=0 = (r-a-l)b0 rf + rf (n+r)("+r+1)+bn+1 (n+r - a)] rf+x. n=0 It requires r=CM-1 and bn+l(n + l) = -bn(n+ a+l)(n + a+2), or (n+a)(n + a~l) b =-b , n n-l n = b0Y[(-m + a/i)(l + (a-l)/i) and (A.7) where b0 =-aa0. The ratio test can be applied to the asymptotic solution (A.5) when a is not a nonpositive integer (otherwise the series is a finite polynomial): lim = Vlim (n + af \ n+l -»0O. 86 The series diverges for any value of r| other than zero, but can provide a good estimate at small values of r|, that is for points in the far-field, under the condition that the series is computed up to the term before the smallest one, say at n=N-l. The latter can be related to the inaccuracy (cf. Murray [1984], section 1.1), and the remaining terms are slowly growing, but they alternate about the estimate. The greater is N, the better is the fit, but the farther the far-field is from the front. The lowest z at which the asymptotic solution can be accurately computed using AT terms is given by \l + (a-l)/N\2\-N/z\ = l. That is (N-1 + Re(a)f +Im(af z = • N The lowest z for which the derivative of the asymptotic solution can be accurately computed up to the JV6 term is given by |1 + (a - l)/N\\l + a/N\\- N/z\ = 1, or z = N\] . 2a-l a(a-l) 1+ + v ' N N2 The appropriate z0 at which the asymptotic expansion of the Tricomi function and its derivative are accurate with at least N terms (this criterion will mark the far-field) is z 0 = y V m a x 1 + a-\ N 1 + a-1 N 1 + If the transformation of the endpoint of the front gives a z = 2K//l that is greater or equal to z 0 (the wavenumber being high enough or the slope small enough to locate the endpoint in the far-field), then the boundary conditions are directly computed by the far-field expansion (A. 5, A. 7). Otherwise, if z is small enough (near-field), the characteristic solutions of the expansion about the singular point are used with the appropriate matching to the far-field boundary conditions. Again, the asymptotic expansion applies to a case where the slope stretches indefinitely. In the case of an horizontal bottom away from the endpoint of the front, z 0 is reset as the 87 location of the foot of the ridge. After this point, the left-hand side of the continuity equation (1.11) is zero. Thus the off-ridge pressure perturbation field is simply exponentially damped: Thus p is a constant beyond the foot of the ridge, say unity, and q=0. The boundary conditions at z=z0 are now p-\, q=0. An overlap region between the far-field and the near-field domains does not necessarily exist for all values of a; z may be neither small enough to be located in the near field, nor large enough to be in the farfield. The boundary conditions are computed by use of analytical continuity, with a chain of matching points relating intermediary expansion solutions about ordinary mid-field points to the asymptotic solution and its derivative (or to the flat bottom solution) at the point z0: where zN is chosen as the endpoint of the front if greater than one (say in the mid-field), or, if the endpoint of the front is in the near-field, zN is chosen as unity and independent solutions of the expansions about z0 are matched to the solutions at zN=l, before their evaluation at the endpoint of the front. For an expansion about the ordinary point z, of the mid-field, let x=z-zt. The ordinary differential equation for p becomes pcce z>z0. (p,q)M^> (p,q)(zi )->•••-> (p,g)(zN) (x + z, )p" + ( l - zt -x)p'-ap = 0, and the one for q becomes (x + zt )q" + (2 - z, - x)q' -(a+ l)q = 0. The series expansions for p and q are anx ,n The substitution of p in its ordinary differential equation leads to 0 = X(* + Z,X"(W-l)x"~2 + £ (l - x - z,)an x"~x - a]Ta n x n=0 n-0 n=0 n=2 n=\ n=0 = X [(" + 2)(« + l)ztan+2 +(n +1)(« +1 - z, )aH+1 - (n + or)a„ ]x The validity for any x in the neighborhood of z, requires the indicial equation: a = (" + « K - (» + +1 - -g/ K + i ( w > 0 ) fl+2 (w + l)(« + 2)z, Similarly for we get 0 = £ [(n + 2){n + l)z,A, + 2 + (it +1)(« + 2 - z, K + 1 - (« +1 + ]: n=0 and its indicial equation is K = q(zi\ a/?(z,.) + (z,-l)tf(z,.) * i = : > 89 APPENDIX B This section contains the subroutines pertaining to the search for eigenvalues in the complex plane. The main file must look like: #include <stdlib.h> #inolude <stdio.h> #include <math.h> void main (), scant); byte bissect2(), recorditO; int isgn() ; void main() { int m, 1; FILE *res; res= fopen("result.txt","a+"); initialise(); //configure the model for(K=6.; K<=9.5; K+=.2) { File. count=0; update () ; scan () ; for(m=0;m<File.count; m++) fprintf(res,"%7.4f,%12.6e,%12.6e,%12.6e\n", K, File.table[m].x, File.table[m].y, File.score[m]); } fclose(res); > // new types are defined typedef unsigned short byte; typedef struct POINT (float x,y;} point; typedef point *p_point; typedef struct SIGNS (int r, i;} signs; typedef signs *p_sign; // global variables byte ALERT; short F; float W, EPS=l.e-7F; point GRID; struct ARCHIVES (int count; double K; point table[MX]; float score[MX];} File; The subroutine scan runs over a given rectangular domain of the complex domain to fill, for a given frequency, an array of signs for the real and imaginary parts of the dispersion function at different values of growth rate. When the array is filled or the lower bound is reached or computational hassles are signaled by the flag ALERT, a new array is filled at a nearby frequency. At every level, signs are compared with the previous level and those of the previous array. When simultanous sign changes occur, a bidimensional bissection technique is applied. Finally the second array is transferred to the first and scanning is reiterated. #define maxi 40 //allowed size of arrays void scan() { float wr, wi, wrd, w=.95F, x, y, ro, border, top, dwr, dwi; int i , j , count; dcomplex z; signs a[maxi], b[maxi]; short int k, 1; ro=(float)Ro;//Rossby # top= 0.01235F; y= 0.0115F;//upper and lower growth rates (example) wr= w*ro; border= -w*ro;//upper and lower frequency bounds ALERT=0; z=dispersion(wr, top); wi=top; b[0].r= isgn(z.r); b[0].i= isgn(z.i); // get signs of dispersion function if(!b[0].rs&!b[0].i) recordit(wr,wi); for (i=l; i<maxi; i++) {//a f i r s t array i s computed wi*= w; z=dispersion(wr, wi);//decrease exponentially the growth rate if(!ALERTCSwi>=y) { b[i].r=isgn(z.r); b[i].i=isgn(z.i); if(!b[i].r&£!b[i].i) recordit(wr,wi);//directly on a root } else {count=i;break;} } if(!ALERT) count=maxi; 90 while (wr>border) {// next array i s computed dwr=GRID.x=fabs(wr)?(fabs(wr)>ro/35 ? fabs(wr)/15 : fabs(wr)) : ro/35; wrd= wr-dwr; ALERT=0; z=dispersion(wrd, top); wi=top; a[0]=b[0]; b[0].r= isgn(z.r); b[0].i= isgn(z.i); if<!(b[0].r||b[0).i)) recordit(wrd,wi); for (i=l; i<count; i++) { GRID.y= dwi=wi/20; wi*=w; z«dispersion(wrd, wi); a[i]=b[i]; i f ( ! ALERTS&wi>=y) { b[i].r==isgn(z.r); b[i].i=isgn(z.i); if(b[i].r||b[i].i) < j - i - l ; lc=abs (b[j] .r+b[i] .r+a[i] .r+a[j] .r) ; l=*abs (b[j] -i+b[i] .i+a[i] .i+a[j] .i) ; if(k<3&&l<3) {//check for simultaneous changes of signs if<(a[j].r||ati).i)S&(a[i].r||a[i].i)&S(b[j].r||b[j].i>) { F=0; bissect2(Sb[i],£b[j],«a[i],Sa[j],fiwrd,Swi,dwr,dwi); } } } else recordit(wrd,wi); / / d i r e c t l y on a root ALERT=0; } else {covint=i;break; } / / i f ALERT occurs } /*end for*/ wi=top*pow(w,count-l)///decrease exponentially the growth rate if(!ALERT) for (i=eount; i<maxi; i++) { wi*=w; z= dispersion(wrd,wi); if(!ALKRT&&wi>y) { b[i].r=isgn(z.r); b[i].i=isgn (2.i); if(!b[ij.rfi&«b[i].i) recordit(wrd,wi); } else { covint=i;break; } } wr= wrd; ) } #undef maxi byte recordit(wr, wi) The subroutine accepts or rejects a new root, compares, and then stores it in a global array. float wr, wi; { double score; int i ; byte found; score= DCabs(dispersion(wr,wi)); i= File.count-1; found= (score<.5) ? 2 : 0;//check i f i t i s a pole while(i>=0 && fovjnd==2) {// compare with previous roots i f ( fabs(File.table[i] .x - wr) <= GRID.x && fabs(File.table[i]-y - wi) <= GRID.y ) { if(score<File.score[i]) { File.table[i].x=wr; File.table[i].y=wi; File.score[i]=score; found= 1; } else found= 0; } — i ; } i f (found=2) {//store i n an array File.table[File.count].x= wr; File.table[File.count].y= wi; File.score[File.count]= score; ++File.count; } i f (File. count—MX) { char c; printf ("Pause: MX\n") ; scanf ("%c" ,c) ; } return (found!=0); > byte bissect2(ul, dl, ur, dr, x, y, dwr, dwi) The subroutine bissect2 is a recursive scheme applying a bissection in two dimensions. The inputs are pointers to the signs found at the upper left, lower left, upper right and lower right corners of the grid, the pointers to given midpoints and grid size lengths. The global variable F indicates if it shall be an horizontal or vertical bissection. p_sign ul, ur, dl, dr; float dwr, dwi, *x, *y; { signs a, b; float w; static dcomplex z; static short int k, 1 ; if(!F) F=l; i f (F===l) { // bissect v e r t i c a l l y dwi/=2; w=*y; w+=dwi; W=*x; W+=dwr; z= dispersion(*x,w); a.r= isgn(z.r); a.i= isgn(z.i); if('ALERT) (if(>(a.r||a.i)) (recordit(*x,w); goto END;}} else goto END; z= dispersion(W,w); b.r= isgn(z.r); b.i= isgn(z.i); if(!ALERT) (if(!(b.r||b.i)) (recordit(W,w); goto END;}} else goto END; // upper part k=abs(ul->r+ur->r+a.r+b.r); l=abs(ul->i+ur->i+a.i+b.i); if(k<4&&l<4) {//Bissect2 i s reapplied on upper h a l f - c e l l , if(dwr>EPS*max(fabs(*x),6RID.x)) {//horizontal bissection. F=-l; bissect2(ul,&a,ur,Sb,x,Sw,dwr,dwi); F=l; } else ( dwi > EPS*max(w,GRID.y))//vertical bissection. ? bissect2(ul,&a,ur,&b,x,£w,dwr,dwi): recordit(*x,w); ) // lower part k=abs(dl->r+dr->r+a.r+b.r); l=abs(dl->i+dr->i+a.i+b.i); if(k<4&&l<4) {//Bissect2 i s reapplied on lower h a l f - c e l l , if(dwr>EPS*max(fabs(*x),GRID.x)) {//horizontal bissection. F=-l; bissect2(&a,dl,Sb,dr,x,y,dwr,dwi); F=l; } else ( dwi > EPS*max(*y,GRID.y))//vertical bissection. ? bissect2(&a,dl,£b,dr,x,y,dwr,dwi): recordit(*x,*y); } } else { // bissect horizontally ( F==-l ) dwr/=2; w=*x; w+=dwr; W=*y; W+=dwi; z- dispersion(w,*y); a.r= isgn(z.r); a.i= isgn(z.i); if{!ALERT) {if(!(a.r||a.i)) {recordit(w,*y); goto END;}} else goto END; z«= dispersion(w,W) ; b.r= isgn(z.r); b.i= isgn(z.i); if(!ALERT) {if(?(b.r||b.i)) (recordit(w,W) ; goto END;}} else goto END; // l e f t part k=abs(ul->r+dl->r+a.r+b.r); l=abs(dl->i+ul->i+a.i+b.i); if(k<4&&l<4) {//Bissect2 i s reapplied on l e f t h a l f - c e l l , if(dwi>EPS*max(*y,GRID.y)) {//vertical bissection. F=l; bissect2(ul,dl,sa,Sb,x,y,dwr,dwi); F=-l; } else { (dwr>EPS*max(fabs(*x),GRID.x))//horizontal bissection. ? bissect2(ul,dl,&a,&b,x,y,dwr,dwi): recordit(*x,*y);} } // right part k=abs(dr->r+ur->r+a.r+b.r); l=abs(ur->i+dr->i+a.i+b.i); if(k<4&&l<4) {//Bissect2 i s reapplied on right h a l f - c e l l , if(dwi>EPS*max(*y,GRID.y)) {//vertical bissection. F=l; bissect2(&a,Sb,ur,dr,Sw,y,dwr,dwi); F=-l; } else { (dwr>EPS*max(fabs(w), GRID.x))//horizontal bissection. ? bissect2(Sa,&b,ur,dr,Sw,y,dwr,dwi): recordit(w,*y);} } } END: if(ALERT) printf("ALERT"); return; } int isgn(x) double x;{ return (x>0.) ? 1 : {(x<0.) ? - 1 : 0 ) ; } 92 This section contains the subroutines pertaining to the calculation of the dispersion function. They use predefined functions for complex arithmetics: double Pyth(a,b); // Euclidean norm of a vector handling unbounded numbers double DCabs(); // same function as cabs(), but i t uses Pyth(). dcomplex dcmplx(); // convert two doubles to a dcomplex number dcomplex DCadd(); // add two dcomplex numbers dcomplex DCsub(); // substract second dcomplex number from f i r s t dcomplex DCmul(); // multiply two dcomplex numbers dcomplex DCdiv(); // divide f i r s t dcomplex number by second dcomplex RDCdiv(); // divide double by dcomplex dcomplex RDCmulO; // multiply double by dcomplex The computing subroutines are simplified by defining structure types and constants: typedef struct DCOMPLEX < double r, i ; } dcomplex; typedef struct PQ {dcomplex p, q;} pq; typedef dcomplex matrix[5][5]; typedef matrix *p_mat; typedef dcomplex eigensol[3][2][3]; typedef eigensol *p_sol; const dcomplex ZERO={0.,0.},ONE={l.,0.}; const pq P_Q,={1. ,0. ,0. ,0. } ; Macros for fast inversion or reflection of dcomplex: dcomplex CXI, CX2, CX3, CX4; double AT, BT; #define SQRNRM(a) (CXl=(a), CXI.r*CX1.r+CXl.i*Cl.i) #define CINV(a)(CX2=(a),AT=SQRNRM(CX2),dcmplx(CX2.r/AT,CX2.i/AT)) #define CNEG(a) <CX4=(a), CX4.r=-CX4.r, CX4.i=-CX4.i, CX4) Global variables and subroutines declarations: double C[2][2], 61, G2, Ri, YY1, YY2, YY3, Ro, K, KK, BB1, BB2, ZZla, ZZlb, ZZ2a, ZZ2b, FACTOR, PI; / / I n i t i a l i s e d with the configuration dcomplex dispersion(), determinant{), Ll, L2, Bl, B2; pqbc(), intermediate() , plateau() , PQ1, PQ.2; void matchings () , series (), swapitO, update (), initialise(); initialise()//setup the configuration parameters { Ri= 66.; 61=.827;G2=1.12;//Richardson number and slopes YY2=5.; YY1=2.333;//distances from front endpoints to ridge feet //YYl={0:shelf-break;>0: foot at y=-YYl-YY3} //YY2={0:false bottom;>0: foot at y=l+YY2} YY3=0; //length of the plateau PI=acos(-l.); BB2=BB1=-Ri; BB2/=1.+G2; //prepare modified Richardson numbers ) update()//readjust the Rossby number, the transformations of the path length over the plateau (FACTOR) and of the positions of the front endpoints and ridge feet over the slopes (ZZla,ZZlb,ZZ2a,ZZ2b)for a new wavenumber K. { double k2; k2= 2*K; KK=K*K; Ro=K/(2*Ri); FACTOR=exp(K*YY3); File.K=K; if(61>0.) { ZZla= k2/61; ZZlb= ZZla + k2*YYl; } if(G2>0.) { ZZ2a= k2*(1.+62)/G2; ZZ2b*= ZZ2a + k2*YY2; } > dcomplex dispersion(wr,wi) This function takes the iterate of the complex frequency and returns a relative value of the objective function associated with the dispersion relation. From the set of configuration parameters of the model, it infers another set of parameters (see section 1.4 of the thesis). With those parameters, it invokes the subroutine BC for the solution of the confluent hypergeometric equations which provide the conditions for both front boundaries (unless it is a flat bottom case). Also, it invokes the subroutine series to compute the eigenfunctions of the frontal equations at half the frontal length. The eigensolutions are returned in the array p(i=l,2; r=l,2; j=0,2), where i=l,2 is the layer index, r=l,2 is the BC-tied/free solution index, and j=0,2 is the 93 order of the derivative. The objective function associated with the dispersion relation (in short: the dispersion function) is the determinant of the matrix of the system (1.53) renormalized by the highest norm from each column and each row, and by the highest element in the last addition. In this way, unbounded terms are avoided and the magnitude measures the closeness of the iterate to an eigenvalue. When computation is not possible (as in the vicinity of a pole), the procedure aborts after modifying the global variable ALERT * from 0 to -1,1,2 or 3. float wr, wi; inputs: components of the complex iterate of the eigenvalue. The real part is the nondimensional Doppler-shifted phase frequency and the imaginary part is the growth rate. { double w, y=.5, rr; dcomplex wl, w2, alpha; pq qp; matrix M; int i , j ; C[0][0]=C[0][1]=C[1][0]=C[1][1]= 1.; //scale factors of eigensolutions ALERT= 0; w= fabs((double)wr); i f ( wi<2e-5 SS (fabs(w/Ro-1) < 4e-2 || fabs(fabs(w-1)/Ro-1) < 6e-4) ) { ALERT=1; return; } // COx,G)2: wl= dcmplx((double)wr+Ro,(double)wi); w2= dcmplx((double)wr-Ro,(double)wi); // px,fi2,Xx,X2: Bl= DCmul(wl,wl); Bl.r -=1.; Bl.r*=BBl; Bl.i *=BB1; B2= DCmul(w2,w2); B2.r -=1.; B2.r*=BB2; B2.i *=BB2; Ll= RDCdiv(-K,wl); L2= RDCdiv(K,w2); Ll.r+=Bl.r; Ll.i+=Bl.i; L2.r+=B2.r; L2.i+=B2.i; qp= P_Q; if(61 != 0.) { //Compute (A>9I) the l e f t boundary conditions. alpha=RDCdiv(-l,wl); alpha.r=++alpha.r/2; alpha.i/= 2; qp= be(alpha, ZZla, ZZlb, qp); if(ALERT=5) printf ("Alert 5 at bc#l for w=(%g, %g) ", wr, wi) ; } qp.q.r= qp.p.r-2*qp.q.r; qp.q.i= qp.p.i-2*qp.q.i; if(YY3 != 0.) qp= plateau(qp); qp.q.r*= K; qp.q.i*= K; rr= Pyth(Pyth(qp.p.r,qp.p.i), Pyth(qp.q.r,qp.q.i)); PQ1.p.r= qp.p.r/rr; PQ1.p.i= qp.p.i/rr; PQl.q.r= qp.q.r/rr; PQl.q.i= qp.q.i/rr; qp= P_Q; if(G2 != 0.) {//Compute {p2,q2) the right boundary conditions. alpha=CINV(w2); alpha.r=++alpha.r/2; alpha.i/«2; qp= be(alpha, ZZ2a, ZZ2b, qp); i f (ALERT=5) printf ("Alert 5 at bc#2 for w=(%g, %g) " , wr, wi) ; } qp.q.r= K*(qp.p.r-2*qp.q.r); qp.q.i= K*(qp.p.i-2*qp.q.i); rr= Pyth(Pyth(qp.p.r,qp.p.i), Pyth(qp.q.r,qp.q.i)); PQ2.p.r=qp.p.r/rr; PQ2.p.i=qp.p.i/rr; PQ2.q.r=qp.q.r/rr; PQ2.q.i=qp.q.i/rr; matchings(y, SM);//compute solutions of the coupled system and establish the matrix of matching conditions at y=l / 2 . for (i=l; i<=4; i++) { for (j=l,rr=0.; j<=4; j++) rr=max(Pyth(M[i][j].r,M[i][j].i),rr); if(rr>0.) for(j=l;j<=4;j++) { M[i][j].r/=rr; M[i][j].i/=rr; } }//rescale rows of the matrix of matching conditions return determinant(SM);//its determinant i s the dispersion function J pq plateau (qp) This function modifies the left boundary conditions for a blunt crest using the global variable FACTORS exp(K*YY3) pq qp; { dcomplex a, b; 94 qp.q.r/= K; a.r=(qp.p.r + qp.q.r)/2; b.r= (qp.p.r - qp.q.r)/2; qp.q.i/= K; a.i=(qp.p.i + qp.q.i)/2; b.i= (qp.p.i - qp.q.i)/2; qp.p.r= a.r*FACTOR +b.r/FACTOR; qp.p.i= a.i*FACTOR 4-b. i /FACTOR; qp.q.r** a.r*FACTOR -b.r/FACTOR; qp.q.i= a.i*FACTOR -b.i/FACTOR; return qp; ) void matchings(y, M) This subroutine computes the series expanded from the left frontal endpoint at y, and those expanded from the right frontal endpoint at 1-y with a swap of indiced parameters. double y; p_nvat M; { eigensol p; int k=0; series(PQ1, k, y, &p); /* B.C.-attached solution on the left*/ (*M)[1][1]= p[l][0][0]; (*M)[2][1]= p[l][0][l]; <*M)[3][1]= p[2][0][0]; (*M)[4][1]= p[2][0][l]; /* free part of solution on the left*/ (*M)[1][2]= p£l][l][0]; (*M)[2][2]= p[l][1][1]; (*M)[3][2]= p[2][l][0]; (*M) [4] [2]= p[2] [1] [1] ; i f (ALERT"=0) { k<=l; swapi t(&Bl,6B2);swapit(SLl,£L2); series(PQ2, k, y, fip); /* B.C.-attached solution on the right*/ (*M)tl][3]= CNEG(p[2][0][0]>; (*M)[2][3]« p[2][0][1]; (*M)[3][3]= CMEG(p[l][0][0]); (*M)[4][3]= p[l][0][1]; /* free part of solution on the right*/ (*M)[1][4] = CNEG(p[2][1][0]); (*M)[2][4]= p[2][l][l]; <*M)[3][4]= CNEG(p[l][1][0]); (*M)[4][4]= p[l][1][1]; swapit(&B1,&B2); swapit(SL1,6L2); } } void swapit(a,b) //swap addresses of indiced parameters dcomplex *a, *b; { dcomplex tmp; tmp= *a; *a=»*b; *b= tmp; } void series(qp, k, y, p) The eigensolutions of the frontal equations are solved by power series, for the global parameters KK, Ll, L2, Bl, B2, the boundary conditions enclosed in qp and at a distance y from the boundary inside the frontal region. The eigensolutions are returned in the array p(i=l:2;r=0:l;j=0:2), where z'=l,2 is the layer index, r=l,2 is the BC-tied/free solution index, and/=0,1,2 is the order of the derivative. pq qp; int k; double y; eigensol *p; { int r; register int i , j ; unsigned int n, m, Maximum=240; eigensol BC; FILE *output; double y2, dl, d2, d3, Eps=2e-7, Dl, D2, greatest, biggest; dcomplex c l , c2, c3, c4, c5, a[4], b[3], u[4], eqn; pq ode; These parameters stay constant inside the subroutine: y2=y*y; dl= -KK*y; d2= KK*y2; d3= l.-y; cl.r= Ll.r+KK; cl.i== L l . i ; c2.r= L2.r*y; c2.i»= L2.i*y; e3.r= -B2.r*y; c3.i= -B2.i*y; c4.r= -Ll.r-KK*d3; c4.i= - L l . i ; c5.r= dl-L2.r; e5.i= -L2.i; BC[1][0]C0]=qp.p; BC[1][0][l]=qp.q; BC[1][0][2]= ZERO; // BC-linked solution BC[1][1][0]=ZERO; BC[1][1][l]=ZERO; BC[1][1][2]= ONE; // free solution The BC-linked solution is constructed first (r=0), then the free solution (r=l). for (r=0; r<2; r++) { 95 ode.p.r= el.r*BC[l][r][0].r-cl.i*BC[l][r][0].i+BC[l][r][1].r-2*r; ode.p.i>= cl.r*BC[l][rj[0].i+cl.i*BC[l][r][0].r+BC[l][r][1].i; The boundary conditions for the second pressure field: BC{2][r][0]= DCdiv(ode.p, Bl); BC[2][r][1].r= L2.r*BC[2][r][0].r - L2.i*BC[2][r][0].i - B2.r*BC[l][r][0].r + B2.i*BC[l][r][0].i; BC[2][r][1].i= L2.r*BC[2][r][0].i + L2.i*BC[2][r][0].r - B2.r*BC[l][r][0].i - B2.i*BC[l][r][0].r; ode.p= DCmul(L2, BC[2] [r] [1]) ; ode.q= DCmul(B2, BC[1][r][1]) ; BC[2][r][2].r= (ode.p.r + KK*BC[2][r][0].r - ode.q.r)/4; BC[2][r][2].i= (ode.p.i + KK*BC[2][r][0].i - ode.q.i)/4; The 3 next terms of the first series and 2 next terms of the second series: a[3].r= C[k][r]*BC[lJ[r][0].r; a[3],i= C[k][r]*BC[l][r][0].i; a[2].r= C[k][r]*BC[l][r][1].r*y; a[2].i= C[k][r]*BC[l][r][1].i*y; a[l].r= C[k][r]*BC[l][r][2].r*y2; a[l].i= C[k][r]*BC[l]tr][2].i*y2; b[2].r= C[k][r]*BC[2][r][1].r*y; b[2].i= C[k][r]*BC[2][r][1].i*y; b[l].r= C[k][r]*BC[2][r][2].r*y2; b[l].i= C[k][r]*BC£2][r][2].i*y2; (*P)tl][r][0]= DCadd( DCadd(a[l],a[2]), a[3]); (*P)[X][r][2]= RDCmul(2./y2, a[l]); <*P)tl]Ir][1].r- (2*a(l].r+ a[2].r)/y; (*P)tl][r][1].i= (2*a[l].i + a[2].i)/y; (*P)12][r][0].r= btl].r+ b[2].r + C[k][r]*BC[2][r][0].r; (*P) [2] [r] [0] .i= b[l].i+ b[2].i + C[k][r]*BC[2][r][0].i; (*p)[2][r][2]= RDCmul(2./y2, b[l]); (*P)[2]tr][l].r= (2*b[l].r+ b[2].r)/y; < * P ) [ 2 ] [ r ] ( 2 * b [ l ] . i + b[2].x)/y; n=2; biggest= 1/Eps; Dl= n/y; while(biggest>Eps && n<Maximum) { //compute next terms u n t i l convergence or Maximum terms reached. //Dl and D2 are f i r s t and second derivative operators on an n t h power. m=n; D2= Dl; Dl= (++n)/y; D2*= Dl; a[0].r= m*a[l].r/Dl + ( d.r*a[2].r - cl.i*a[2].i - Bl.r*b[2].r + Bl.i*b[2].i + dl*a[3].r )/D2; b[0].r= ((o2.r*b[l].r + d2*b[2].r + c3.r*a[l].r - c2.i*b[l].i - c3.i*a[l].i)/n)/n; a[0].i= m*a[l].i/Dl + ( cl.r*a[2].i + cl.i*a[2].r - Bl.r*bt2].i - Bl.i*b[2].r + dl*a[3].i )/D2; b[0].i=* ((o2.r*b[l].i + d2*b[2].i + e3.r*a[l].i+c2.i*b[l].r + c3.i*a[l].r)/n)/n; (*P) [1] [r] [0] .r +- a[0] .r; (*p) [1] [r] [0] . i += a[0] . i ; (*P) II]tr][1]-r += a[0].r*Dl; (*p)tl][r][1].i += aJO].i*Dl; (*P) tl]tr][2].r += a[0].r*D2; <*p)[1][r][2].i += afO].x*D2; (*P)12][r][0].r += b[0].r; (*p)[2][r][0].i += b[0].i; (*P)12]tr]tl]-r += b[0].r*Dl; (*p)[2]tr][1].i += btO].i*Dl; (*P)[21[rj[2].r += b[0].r*D2; (*p)[2][r][2].i += b[0].i*D2; // s h i f t previous terms a[3]= a[2]; a[2]=a[l]; a[l]=>a[0]; b[2]=b[l]; b[l]=b[0]; // convergence of series solution are tested i n coupled ODEs. u[0].r= d3*(*p)[1][r][2].r; u[l].r= c4.r*(*p)[1][r][0].r - c4.iM*p)[l][r][0].i; u[2].r= -(*p)[1][r][1].r; u[3].r= Bl.r*(*p)[2][r][0].r - Bl.i*(*p)[2][r][0].x; u[0].i= d3*(*p)[1][r][2].i; u[l].i= c4.r*(*p)[1][r][0].i + c4.i*(*p) [1] [rj [0] .r; u[2].i= -<*p)[1][r][1].i; u[3].i= Bl.r*(*p)[2][r][0].i + Bl.i*<*p)[2][r][0].r; for(eqn=ZERO, ode.p=ZERO, i=0; i<4; i++) { eqn.r= max(eqn.r, fabs(u[i].r)); ode.p.r += u[i].r; eqn.i= max(eqn.i, fabs(u[i].i)); ode.p.i +=u[i].i; } greatest= Pyth(eqn.r,eqn.i); ode.p.r/= eqn.r>0. ? eqn.r : 1; ode.p.i/= eqn.i>0. ? eqn.i : 1; biggest= max(fabs(ode.p.r),fabs(ode.p.i)); u[0].r= (*p)[2][r][2].r*y; u[0].i= (*p)[2][r][2].i*y; u[l].r= c5.r*(*p)[2][r][0].r - c5.i*(*p)[2][r][0]-i; 96 u[l].i=* c5.r*(*p)[2][r][0].i + c5.i*(*p)[2][r][0].r; u[2].r= (*p)[2]tr][1].r; u[2].i= <*p)[2][r][1].i; u[3].r= B2.r*(*p)[1][r][0].r - B2.i*(*p)[1][r]t0].i; u[3].i= B2.r*(*p)[1][r][0].i + B2.i*(*p)[1][r][0].r; for(eqn=ZERO, ode.q=ZERO, i=0; i<4; i++) { eqn.r= max(eqn.r, fabs(ufi].r)); ode.q.r += u[i].r; eqn.i= max(eqn.i, fabs(u[i].i)); ode.q.i += u[i].i; } greatest^ max(greatest, Pyth(eqn.r,eqn.i)); ode.q.r/= eqn.r !=0. ? eqn.r : 1; ode.q.i/== eqn.i !=0. ? eqn.i : 1; biggest= max(biggest, max(fabs(ode.q.r), fabs(ode.q.i)) ); if(max(l.,D2)*DCabs(a[0])<le-15*min(DCabs((*p)[X][r]t0]), min(DCabs((*p)tl][r]tl]),DCabs((*p)[1][r][2]))) SS max(l.,D2)*DCabs(b[0])<le-15*min(DCabs((*p)[2]tr][0]), min(DCabs((*p)t2][r][1]),DCabs((*p)12][r][2])))) break; if(greatest>le9) if(y=.5) { If any solution becomes very large, then we renormalize each accumulated term by greatest to avoid the numeric ceiling. Reseating must be including for any profiling. Ctk]tr] /= greatest; for(i=l; i<=2; i++) for(j=0; j<=2; j++) { (*p)ti]tr]tj]-r/=greatest; (*p)ti][r][j].i/=greatest; } for(j=0; j<=3; j++) {a[j].r /= greatest; a[j]-.i /= greatest;} for(j=0; j<=2; j++) {b[j].r /= greatest; b[j].i /= greatest;} } > } if(n>Maximum){ALERT=2; printf("Alert:2 too many terms needed");} dcomplex determinant(M) This subroutine tests solutions for stiffness, computes the determinant and normalises according to its greatest entry before the last adding for pole detection. matrix M; { double rr, i i , r i , distinct; dcomplex det={100.,100.},ut6],df6],v[6]; register int j ; FILE *output; of the v[0]x v[l] = v[2] = v[3] = vt4] = vt5] = for(j=0, distinct=ri=0.; j < 6; j++) { rr=max(fabs(u[j].r),fabs(v[j].r)); ii=max(fabs(u[j].i),fabs(v[j].i)); ri= max(ri, max(Pyth(u[j].r,u[j].i),Pyth(vtj].r,v[j].i)) ); u[ j] .r+= vt j] .r; u[ j] .i+=v£j] . i ; if(rr=0.) rr= 1.; if(ii==0.) ii= distinct= max(distinct, Pyth(uIj].r/rr, u [ j ] . i / i i ) ) ; } if(distinct<le-6) { // left-hand solutions are s t i f f printf("Alert: 3 at determinant (left)"); ALERT9 3; goto END; } else if(ri>0.) for(j=0;j<6;j++) { u[j].r/= r i ; u[j].i/= ri;} //The f i r s t stage of computation of the determinant M for right series. //The f i r s t stage of computation u[0]= DCmul(M[l][1], Mt2]f2]); u[l]= DCmul(M[3][1], ut2]= DCmul(M[l]tl], u[3]= DCmul(M[2]II], u[4]= DCmul(Mt4][1], u[5]= DCmul(M[3]tl], Mil]t2]) Mt4][2]) M[3]t2]) M[2][2]) Mf4]t2]) determinant M for l e f t series DCmul(CNEG(MI2]tl]), Mil][2]) DCmul(CNEG(Mil][1]), M[3]I2]) DCmul(CNEG(MI4][1]), M[l]t2]) DCmul(CNEG(Mt3]II]), M[2] 12]) DCmul(CNEG(M[2][1]), M[4]I2]) DCmul(CNEG(M[4]tl]), M[3]t2]) dt0]= DCmul(M[3][3] d[l]= DCmul(MC2][3] d[2]= DCmul(Mt2][3] dt3]= DCmul(M[l]13] d[4]= DCmul(M[l]f3] d[5]= DCmul(M[l]13] M[4]t4]); v[0]= DCmul(CNEG(M[4]13]), M[3][4]) Mt4][4]); vtl]= DCmul(CNEG(M[4][3]), M[2]t4]) Mf3][4]); vf2]= DCmul(CNEG(M[3][3]), M[2]t4]) M[4][4]); v[3]= DCmul(CNEG(M[4][3]), M[l][4]) Mt3][4]); v[4]= DCmul(CNEG(M[3][3]), M[l][4]) M[2][4]); v[5]= DCmul(CNEG(M[2][3]), Mtl]f4]) for(j=0, distinct=ri=0.; j < 6; j++) ( rr=max(fabs(d[j].r),fabs(vtj].r)); ii=max(fabs(dtj].i),fabs(v[j] ri=max(ri, max(Pyth(d[j].r,dtj]-i),Pyth(vtj].r,v[j].i)) ); dtj].r+= vtj].r; d[j].i+= v[j].i; if(rr=0.) r r * l . ; if(ii==0.) i i = l . ; i)) 97 distinct= max(distinct, Pyth(d[j].r/rr,d[j].i/ii)); > if(distinct<le-6) {// right-hand solutions are s t i f f printf("Alert: 3 at determinant (right)"); ALERT= 3; goto END; } else if(riX).) for(j=0;j<6;j++) { d[j].r/= r i ; d[j].i/= ri;} //'the second stage involves products of minors for(j=0; j<6;j++) v[j]= DOnul(u[;j] ,d[j]) ; // l a s t stage: the determinant i s summed then normalised det=> ZERO; ii= rr= 0.; for(j=0; j<6;j++) { det.r+= v[j].r; rr= max(rr,fabs(v[j].r)); } for(j=0; j<6;j++) { det.i+= v[j].i; ii= max(ii,fabs(v[j].i)); } if(rr>0.) det.r/=rr; else det.r= 0.; if(ii>0.) det.i/=ii; else det.i= 0.; END: return det; } pg be(alpha, z, zz, qp) Compute the pressure field and its derivative (p,q) over the slope beyond the front through a chain of expansions and matchings from the ridge foot (zz) to the front endpoint (z). dcomplex alpha; double z, zz; pq qp; { double y, x; y= 1.+max(l.,fabs(alpha.r)); //the expansion stepsize i s adapted to alpha x= -rain ( 2 2 - 2 , zz/y) ; while (zz > z) {//chain of intermediate expansions and matchings u n t i l the front endpoint, zz i s the renewed expansion point. qp= intermediate(alpha, fizz, Sx, qp); x=-min(zz-z,2*fabs{x)); //next expansion t r y a greater stepsize } return cjp; > pq intermediate(alpha, zz, xx, qp) This subroutine computes the expansions of the confluent hypergeometric solutions through the midfield. The expansion origin is addressed by zz and the (negative) step is addressed by xx. qp holds the boundary conditions at the expansion origin and will receive the expanded solutions at the next point. The relative accuracy is tested in the ODE. If it is not sufficient, then the separation is halfed and the process restarts with a modified stepsize. dcomplex alpha; double *zz, *xx; pq qp; { double xz, xxz, xzn, xxzn, x, z, z i , z i l , fa, fb, fc, f, g, h, eqr, eqi, norm, eqmx, Eps=2e-7, Inf=le50; dcomplex nalpha, nlalpha, n2alpha, p, q, r, a[3], b[3], c[3], u, v, w, eqn; register int n2, nl2; z- *zz; x=s *xx; BEGIN: zi= z+x; xz= x/z; xxz= x*xz; zil= l . - z i ; p= qp.p; q= qp.q; u= DCmul(alpha, p); r.r= (u.r + (z-1.)*q.r)/z; r.i= (u.i + (z-1.)*q.i)/z; nalpha= alpha; nlalpha.r= l.+alpha.r; n2alpha.r= 2.+alpha.r; nlalpha.i= n2alpha.i= alpha.i; u= DCmul(nlalpha,q); a[0]=p; a[l].r=x*q.r; a[l].i=x*q.i; b[0]=q; b[l].r=x*r.r; b{l].i=x*r.i; c[0]=r; c[l].r=xz*(u.r+(z-2.)*r.r); c[l].i=xz*(u.i+(z-2.)*r.i); p.r+=a[l].r; q.r+= b[l].r; r.r+=c[l].r; p.i+= a[l ] . i ; q.i+= b[ l ] . i ; r.i+= c [ l ] . i ; nl2= n2= 2; fa=l.-z; fb=2.-z; fc=3.-z; eqmx= 1/Eps; while( eqmx > Eps && n2 < 30 ) { xzn= xz/n2; xxzn= xxz/nl2; u.r= a[0].r*nalpha.r - a[0].i*nalpha.i; u.i= a{0].i*nalpha.r + a[0].r*nalpha.i; v.r= b[0].r*nlalpha.r - b[0].i*nlalpha.i; v.i= b[0].i*nlalpha.r + b[0].r*nlalpha.i; w.r= c[0].r*n2alpha.r - c[0].i*n2alpha.i; w.i= c[0].i*n2alpha.r + o[0].r*n2alpha.i; f*=fa*xzn; g=fb*xzn; h=fc*xzn; a[2].r= xxzn*u.r - f*a[l].r; a[2].i= xxzn*u.i - f*a[l].i; b[2].r= xxzn*v.r - g*b[l].r; b[2].i= xxzn*v.i - g*b[l].i; c[2].r= xxzn*w.r - h*c[l].r; o[2].i= xxzn*w.i - h*c[l].i; p.r+= a[2].r; q.r+= b[2].r; r.r+= c[2].r; p.i+= a[2].i; q.i+=b[2].i; r.i+= c[2].i; a[0]=a[l]; a[l]=a[2]; b[0]=b[l]; b[l]=b[2]; c[0]=c[l]; c[l]=c[2]; nl2= r»2++; nl2 *<* n2; nalpha.r= nlalpha.r; nlalpha.r= n2alpha.r; n2alpha.r= alpha.r + n2 fa= fb; fb= fc; fc= (n2+l)-z; u.r= zi*r.r; u.i= z i * r . i ; v.r= zil*q.r; v.i= z i l * q . i ; w.r= alpha.r*p.r-alpha.i *p.i; w.i= alpha.r*p.i+alpha.i *p.r; eqn.r= u.r + v.r - w.r; eqn.i= u.i + v.i - w.i; eqr= max(max(fabs(u.r), fabs(v.r)), fabs(w.r)); if(eqr!=0.) eqr= fabs(eqn.r)/eqr; eqi= max(max(fabs(u.i), fabs(v.i)), fabs(w.i)); xf(eqi!=0.) eqi= fabs(eqn.i)/eqi; eqmx= max(eqr, eqi); norm= DCabs(Pyth(p.r,p.i), Pyth(q.r,q.i)); if(norm>=Inf) break; } if(n2<30 &£ nornKInf) goto END; // i t has not converged yet or the solutions have grown too big x/=2;// i t w i l l be ret r i e d with smaller stepsize if(norm>«=Inf) { qp.p.r/=Inf; qp.p.i/=Inf; qp.q.r/=Inf; qp.q.i/=Inf; } goto BEGIN; END: °P-P'=,P; qp.q= q; *zz= z i ; *xx=x; return qp; > 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items