AN ANALYTIC MODEL OF OBSERVED SMALL SCALE BAROCLINIC CURRENTS IN THE ARCTIC OCEAN Richard Ivan Itkin ^POSTGRADUATE SCHOOL 'iioHTERH. CALIF. 93940 ■:} Monterey, California :. T utm ipmiwm— mini ~ mini AN ANALYTIC MODEL OF OBSERVED SMALL SCALE BAROCLINIC CURRFNTS IN THE ARCTIC OCEAN by Richard Ivan Itkin Thesis Advisor: J. A. Gait September 19/3 Appiovzd (Jo-t pubVLc kzJLq&!>&; cU-A&iibivtion imJUmAnd. An Analytic Model of Observed Small Scale Baroclinic Currents in the Arctic Ocean by Richard Ivan Itkin Lieutenant Commander, United States Navy B.M.E., Cornell University, 1963 Submitted in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE IN OCFANOGPAPHY from the NAVAL POSTGRADUATE SCHOOL September 19 73 • LIBRARY NAVAL POSTGRADUATE SCHOOD MONTEREY, CALIF. 93940 ABSTRACT An analytic model is developed of baroclinic currents measured by several investigators under the Arctic Ocean ice cover. A time-dependent, Boussinesq, quasi-geostrophic vorticity equation is developed to describe the flow. Separable normal mode solutions are found which show sta- tionary current patterns qualitatively consistent with observed data. An initial value problem is then formulated using density perturbations at a newly forming lead as a forcing function. TABLF OF CONTENTS I. INTRODUCTION 8 A. OBJECTIVES OF RESEARCH AND REVIEW OF PRFVIOUS RESEARCH 8 B. THEORY DEVELOPED IN THIS THESIS 20 II. DEVELOPMENT OF THE WAVE SOLUTION ; 2 4 A. ANALYTIC DEVELOPMENT OF THE VORTICITY EQUATION 24 B. EQUATION USING ARCTIC DENSITY PROFILE 31 1. Basic State Description 31 2. Use of /Actual Density Distribution in Equation 33 III. RESULTS OF WAVE SOLUTION AND COMPARISON WITH DATA 37 A. NUMERICAL SOLUTION OF THE VERTICAL MODF EQUATION 37 B. RESULTS OF THEORETICAL DEVELOPMENT 39 C. COMPARISON WITH OBSERVFD DATA : 44 IV. DEVELOPMENT OF THE INITIAL VALUE PROBLEM 51 A. GENERAL METHOD 51 B. SPECIFICATION OF INITIAL CONDITIONS 52 C. SOLUTION OF THE PROBLEM 55 V. CONCLUSIONS 62 VI. RECOMMENDATIONS FOR FURTHFR RFSFARCH 65 BIBLIOGRAPHY 66 INITIAL DISTRIBUTION LIST 68 FORM DD 3 473 70 LIST OF TABLFS Table I. -1 9 Values of a (sec x 10" ) and T (yrs) for mn mn VJ ' Selected Vertical and Horizontal Modes 60 II. Values of k and B for Horizontal Modes m m 1 Through 10 61 III. Values of A and A for Vertical Modes n n 1 Through 6 61 LIST OF DRAWINGS Figure 1. Schematic of Convection Model Proposed by Coachman [1966] 12 2. Proposed Shear Mechanism for Lead Formation 14 3. Vertical Distribution of Horizontal Velocities Measured by Gait [1967] 16 4. Vertical Distribution of Horizontal Velocities Measured by Coachman and Newton [1973] 17 5. Horizontal Distribution of Horizontal Velocities Measured by Coachman and Newton [1973] 18 6. Vertical Section of Velocities Measured by Coachman and Newton [1Q73] 19 7. Observations of Salinity, Temperature, and Density Taken Under Ice Island T-3 22 8. Depth Distribution of Analytic Density Function Showing Relation to Data 32 9. Vertical Structure of Theoretical Horizontal Velocity Patterns for Modes 1 through 4 40 10. Vertical Structure of Theoretical Vertical Velocity Patterns for Modes 1 through 4 41 11. Effect of Layer Depth Change on Vertical Structure 42 12. Effect of Varying e-fold Depth in Density Distribution on Vertical Structure 43 13. Horizontal Distribution of Theoretical Horizontal Velocities 45 14. Linear Sum of Modes 1 and 2 for h=30, c=120 47 15. Vertical Cross-section of Horizontal Velocities, Mode 2 49 16. Vertical Cross-section of Horizontal Velocities, Sum of Modes 1 and 2 50 17. Perturbation Density Due to Salt Rejection at Time Zero 53 18. Vertical Cross-section of Perturbation Density at Time Zero 53 19. Perturbation Pressure at Time Zero ■■ 54 20. Fourier Series Expansion of Initial Horizontal Pressure Distribution 57 21. Fourier Series Expansion of Horizontal Velocity Distribution in x-direction 57 22. Series Expansion of Vertical Structure Function ■ 59 ACKNOWLEDGEMENT The author wishes to sincerely thank Dr. Jerry A. Gait for the many hours spent discussing the formulation and development of the analysis used in this thesis, as well as for the benefit of Dr. Gait's background in the field of Arctic Ocean dynamics. The time spent in careful reading of the manuscript by Dr. R. G. Paquette was also greatly appreciated. This thesis was done in conjunction with a scientific program under the general direction of J. A. Gait, sponsored by the Arctic Section of the Office of Naval Research (Research No. RR132-05-01 NR 307-345). I. INTRODUCTION A. OBJECTIVES OF RESEARCH AND REVIEW OF PREVIOUS RFSFARCH The objective of this study was to develop an analytic model of specific small scale baroclinic current patterns which have been observed on a recurring basis beneath the Arctic ice cap. The patterns are termed small in that they are on a scale less than that of the general Arctic circula- tion. In order to develop the model a theoretical investi- gation was made of the baroclinic response of the Arctic Ocean to perturbations in its internal density structure at a newly forming lead. It is worthwhile at this point to briefly discuss pre- vious research in the area of Arctic dynamic modeling, and to describe the baroclinic currents measured by several investigators, in an effort to indicate the relationship of this study to previous work. Investigations of Arctic Ocean dynamics can be classified by their scale as well as by the mode of motion described (either internal, surface, or a combined mode). The largest scale studied concerns processes which affect the entire ocean, such as long term mean ice drift and water circulation In order to practically study motion on this scale, spacing between grid points cannot feasibly be less than tens of kilometers. This necessitates parameterizing some of the smaller scale processes occurring between grid points, rather than solving for them explicitly. Factors which must be 8 considered on this scale are the shape and bathymetry of the ocean basin, Coriolis effects, wind and ice stress, etc. Several theoretical investigations of large scale Arctic Ocean dynamics which closely match observed patterns of water circulation and ice drift have been carried out in recent years. Among these, Campbell [1965] proposed a large scale analytic ice drift model and Gait [ 197 3a] developed a working barotropic numerical model of overall surface water circula- tion. Both these models were concerned with surface motion or an equivalent barotropic circulation. A complete model of large scale internal /arctic circulation has not yet been published, but work along these lines is being carried out at NOAA's Geophysical Fluid Dynamics Laboratory, Princeton, New Jersey [Semtner, 1973]. The smallest scale of analysis concerns the short term local effects of phenomena such as refreezing leads and ice cracking. These effects result in dynamic processes such as forced convection due to the salt rejection and heat transfer during freezing. Grid size of concern here would be from less than a cm, a microscale, to a few meters. Coriolis effects would not be important at the smallest end of this scale but effects such as diffusion would have to be con- sidered. Several investigators have done work on the problem of forced convection associated with freezing of sea ice. Foster [1972] treated the problem analytically after having published several earlier papers on the problem of forced convection in the ocean. The initial state in his study was a homogeneous under-ice layer. He then assumed density to be a function of salinity only and solved the equations of motion and salt diffusion directly associated with a refreezing lead. Using a two-dimensional model Foster carried his solution to the point of finding initial vertical velocities and length and depth scales for convection cells during the initial stages of convection. These were about 1.4 cm/sec for vertical velocities, 80 cm for horizontal convection cell spacing, and several cms for depth of con- vective penetration. He did not analytically consider further stages of convection nor did he continue to a larger scale in which Coriolis forces would become important. He did postulate a higher order of convection cells with depth penetration of up to 6 meters and vertical velocities of about 1 cm/sec. Piascek [1970] considered the general problem of surface cooling effects on the ocean. He developed a numerical model for convection cells on both a small scale (order of cms) which neglected Coriolis effects and a larger scale (order of kms) which considered Coriolis. His model con- figuration was a rectangular box directly below the surface with an initially homogeneous structure. He considered uniform surface cooling as a boundary condition and then developed a vorticity equation in terms of a vector poten- tial. This was numerically integrated over a 2-dimensional grid to give steady state temperature and vorticity patterns within his "box." With some modification his model could be adapted to describe Arctic under-ice convection. 10 Coachman [1966] considered a mass and energy balance on the layer of water immediately below the ice cover. He assumed a downward flow direction under the refreezing leads due to sinking of highly saline water and a compensating upward flow spread out over the area between leads, as shown in Figure 1. His study used mean values of freezing rate, temperatures, lead size and distance • between leads to deter- mine mean vertical velocity values. His scale of analysis was necessarily larger than that of Foster in order to arrive at mean velocities. His downward vertical velocity values are on the order of .01 to .08 cm/sec. The depth of his surface circulation layer was 100 meters. If the Arctic Ocean is viewed from a large scale standpoint, the net result of the investigations by Coachman, Foster, and Piascek would be to show that there is a small downward vertical velocity associated with each refreezing lead due to forced convection and that there is a local increase in the density of the water directly under a refreezing lead. A related small scale problem which has been investigated is that of the ice formation rate and rates of salt rejection and heat transfer at a refreezing Arctic lead. Schaus [1971] developed a numerical model of a refreezing lead which treats these problems and which could be used in coupling the air- sea interaction problem with ice dynamics and with the under- ice convection problem discussed here. A vertical velocity similar to that caused by the convec- tive process may also arise from the initial cracking and ice motion during lead formation. This specific problem 11 / L !: A D MCE rnnzr- c y w V t / / k j \ v 1 J i \ \ 1 i t Q Figure 1. Schematic of Convection Model Proposed by Coachman [19 6 6] . 12 has not been previously studied. In this case the theory is as follows, assuming ice shears upon cracking over a period of a few days as shown in Figure 2. There is a stress imparted to the underlying water at the ice-water boundary, which in turn causes a convergence of the Fkrnan layer. This in turn results in an impulsive vertical velocity under the newly forming lead. The net result is that, for an entirely different reason, a vertical velocity similar to that caused by convection is associated with shear across newly forming leads. Both these lead-associated phenomena would be classed as small scale effects. By interacting with the existing density and flow patterns these velocity and density perturbations may have some, as yet undetermined, effect on the large scale internal flow dynamics of the Arctic. A proposed theory for this is discussed and developed in part B of this section. Any resultant flow pattern due to such an interaction might then be expected on a scale intermediate in size between the smallest and largest scales discussed; on the order of meters in depth and meters to kilometers in length. A theo- retical investigation of baroclinic internal flow in the Arctic on this scale has not previously appeared in the literature. It is important to note here that this is the scale of a number of internal current patterns or events measured in the Arctic. The observed currents appear to be baroclinic in nature although a satisfactory theory for their origin has not yet 13 ICE MOTION xnzxzszs: WATER MOTION Figure 2. Proposed Shear Mechanism for Lead Formation. 14 been developed. An example of the vertical distribution of horizontal velocities measured by Gait [1967] under Ice Island T-3 is shown in Figure 3. The corresponding density distribution is also shown. Similar results found by Coachman, and Newton [1973] are shovTn in Figure 4. It can be seen that a peak in the horizontal velocity occurs in the pycnocline, with a maximum at about 150 meters, indicating a baroclinic origin. Figure 5 shows the horizontal distri- bution of velocities found by Coachman and Newton at 150 meters. It is a plot of velocity vectors measured from one position on the ice while drifting over a given area. As can be seen, there is a periodic change in the velocity direction with an apparent wavelength of about 30 km. The reason for classing this periodicity as spatial rather than temporal is that the observed currents were found to be essentially in geostrophic balance. For this type feature, phase propagation speeds will be small compared to advective speeds. Figure 6 shows a vertical cross-section of veloci- ties as reported by Coachman and Newton [1973], A need thus exists for an investigation of the Arctic Ocean internal flow dynamics on the intermediate scale described. The immediate purposes of this study are to offer an explanation of the observed baroclinic currents. A further use of such a study is to couple the small scale phenomena to the overall circulation, which in turn may aid in thfi development of a three dimensional Arctic circulation model . 15 100 i- E a a 20 0 - 300 — Velocity (cni/ sec) 4 0 60 80 Figure 3. Vertical Distribution of Horizontal Velocities Measured bv Gait [1967]. 16 Velocity (cm/sec) 10 20 T E a a Figure 4. Vertical Distribution of Horizontal Velocities Measured by Coachman and Newton [19 7 3] . 17 1 w •H -J-) ■H u o rH 0) t-\ en rd r*- +J en o >— ' N •H *H O m o c o T3 fd -P C 3 rd •h jc; ^ u 4-) rd w O •H U C >. rd -P 15 c o; o u N 3 O in Cn ■H P4 IB »■ | 35 Horixonfal Disfance Ckm) 30 25 20 15 T r I 15 1 cm/sec 50 100 ? •u 3" 3 I50 ^ 200 250 300 Figure 6. Vertical Section of Velocities Measured by Coachman and Newton [19 73] . 19 B. THEORY DEVELOPED IN THIS THESIS A theory was postulated to describe the dynamic effects on the underlying ocean of a nev/ly forming Arctic lead. The theory developed was begun with the premise that a small vertical velocity perturbation or density increase is created during the formation of each lead, caused by either refreez- ing or cracking, as described earlier. The scale of analysis used v/as tens of kms in the horizontal and the full depth of water (0-3000 m) in the vertical. The water directly below the Arctic ice cover can typi- cally be expected to consist of a shallow homogeneous layer (20-50 m) , under which a stably stratified density distribu- tion exists. There is also some horizontal slope to the isopycnals in order to support a mean flow. The horizontal gradient will be several orders of magnitude less than that of the vertical density gradient in the pycnocline. Pertur- bations such as the vertical velocity and density changes occurring under a lead, which is a very small area on the assumed length scale, should distort the isopycnals over a relatively short period of time and short distance. This would act as an impulsive event on an overall time scale measured in days. As the isopycnals attempt to return to their stable condition, a transient waveform of some sort would then appear and propagate out from the lead. The first problem would then be to find the form of the waves which can exist in a typical Arctic density structure and to then determine if these are consistent with observed 20 current patterns. Once this was done, a second problem would be an investigation to determine whether the perturbations at a lead could actually force these oscillations. An inves- tigation of these two problems forms the basis of this thesis. An indication of the possible validity of the work at any point is the degree of correlation of theoretically developed waves or currents with real data. A comparison was therefore made with the baroclinic currents measured by Gait [1967] and Coachman and Newton [1973] once a theoretical current pattern was developed. The first step in substantiating the theory was to inves- tigate the type of waves which could exist in the given density structure and to compare these with known data. This corresponded to treating the Arctic Ocean as a semi-enclosed system, subject to certain governing equations of motion and boundary conditions , and solving for its preferred modes of oscillation. The actual solution to this problem is the subject of Section II of this thesis. The general assumptions made in formulating the theory and the rationale behind them will be discussed here. The equations used follow the development of Gait [1973b] which is based on that of Kuo [1973], The basic difference from Kuo's development is in the formulation of the equation of state. The system was assumed to be inviscid, hydrostatic, and quasi-geostrophic . The equation of state assumes density proportional to salinity only. This is a reasonable approx- imation in the region being investigated as can be seen in Figure 7, which is a plot of salinity, temperature, and 21 2800 Figure 7. Observations of Salinity, Temperature and Density Taken Under Ice Island T-3. 2 2 density profiles taken under T-3 along with current measure- ments [Gait, 1967] . The Boussinesq approximation was also made, which states that density variations are significant only when the density is multiplied by gravity. An additional assumption necessary to this development is that the Rossby number is small. This allowed the expansion of the dependent variables as power series in terms of the Rossby number. Equations consistent in powers of Rossby number could then be derived. From the above assumptions and the basic equations of motion a quasi-geostrophic vorticity equation was developed in terms of the perturbation pressure (pressure could be written as the sum of a mean and a perturbation quantity) . This equation was solved by assuming a separable solution consisting of the product of a horizontal plane wave and an unknown vertical structure function. This vertical structure function was then solved for numerically subject to appro- priate boundary conditions. Results of the numerical solu- tion are compared with the measured currents in Section III of this thesis. Investigation of the second problem described (whether the theoretical oscillations could actually be driven by perturbations found at a lead) was accomplished by using pressure conditions under a newly forming lead as an initial condition for the quasi-geostrophic vorticity equation developed in Section II. A description of this initial value problem forms Section IV of this thesis. 23 II. DEVELOPMFNT OF THE WAVE SOLUTION A. ANALYTIC DEVELOPMENT OF THE VORTICITY EQUATION The initial formulation of the problem was carried out using a rectangular coordinate system, with positive values of x to the East, y to the North, and z upward. The govern- ing equations of motion for the flow can be written in the following Boussinesq, inviscid form: u. + uu + vu + wu - f v = - — p (1) t x y z p x x 1 s v, + uv + w + wv + f u = - — p (2) t x y z ps y w, + uw + vw + ww = - — p - — g (3) t x y z p *z p ^ 2 s s where u,v, and w are the x,y, and z velocity components, respectively, f is the vertical component of the Coriolis parameter, g is the acceleration of gravity, p is the pres- sure, p is the density and p is the horizontal space and time average of the density, which is a function of z only. In addition to the equations of motion, the continuity equa^ tion can be written as: p.+up +vp + wp +p(u+v+w)=0 (4) t x y z s x y z and the equation for conservation of salt, assuming a con- servative lav; with no mixing, is: s. + us + vs + ws = 0 (5) t x y z 24 Equations (1) through (5) can then be scaled and non- dimensionalized using essentially the same system as was used by Kuo [1973]. The non-dimensional quantities (primed) are defined as follows: (x,y) = L(x',y'); (u,v) = U(u',v'); (6) z = Dz ; w = — - — w ; \ L I t = (16) t x y z ^x R(v^ + uv + vv + Rwv ) + fu = -4 (17) t x y z y + 6 6 + s = 0 (18) Yz pY u + v + — (wp ) =0 (19) x y ps s'z R(s. + us + vs + Rv;s + -2- w) = 0 (20) t x y z 1 \i I At this point the assumption of a small Rossby number, R, will be used and dependent variables expanded in a power series form; i.e., u, v, 6, and s will take the form: a = E Rna (21) n=0 n where the a represents the general dependent variable and w will be written as: 2 w = w, + Rw„ + R'w-, + ... etc. This is consistent with the scaling given by equation (6). In addition, the Coriolis parameter will be scaled on a beta-plane : fQ = 1 + By (22) where By < 0(R) (23) Using these approximations plus eg. (15) and the series 27 expansion gives the following equations to zero order in Rossby number: -v = -(J) (24) o o X u = -4> (25) o Yo Y = -s (26) o o z u + v =0 (27) o o x y These show the flow to be geostrophic, hydrostatic, and hori- zontally non-divergent to a zeroth order approximation. The equations to first order in Rossby number are: u +uu +vu - v, - 3yv = -d), (28) o, o o o o 1 J o Y] t x y x v + u v +vv + u, + 3yu = -, (29) o , o o o o 1 2 o ^ ly t x y J ♦l ■ "Vo - s.l (30) z K ui + vi + r + u ■? — 6 + v -? — - o o , o I o o o 1 6 o , (36 p z' t P z X P z' y and differentiating this with respect to z gives: (wlpsl = /yps \ fyps 1 7 d) + U 7 d) + V 0 o / . o o o / o P z' zt l P z' zx yp. 6 Yo p z zy + Uo 6—^0 + Vo z ' p z' x : /ypQ o o p .z' (37) y which can then be substituted directly into (34). In writing this it is helpful to define all of the velocities in terms of the pressure, assuming a geostrophic relationship (this assumption filters out gravity waves) . u = -d> o ro y (v -u ) = d> + d) °x °y XX ^ o Yo X (38) From these it can be seen that the last two terms in (37) 29 identically cancel. Then, substitution of (37) and (38) into (34) and dropping all zero subscripts gives: ^t x^y q = 0 Yynx (39) where : 'yp, q = A + (J) + By - — M vxx Yyy M^ ps (40) Equation (39) is the quasi-geostrophic vorticity equa- tion and (40) is the potential vorticity. The first two terms on the right hand side of (40) constitute the relative vorticity. The third term is the contribution due to the variation in the Coriolis parameter, and the last term repre- sents the vorticity contributed by vertical stretching of the water column. Equation (39) shows that time-dependent motion will result when the motion is such as to create a change in potential vorticity. For the initial problem studied here, the Coriolis parameter, 6, was assumed to be non-zero and the resultant motion investigated. It should be noted that as the North pole is approached, 3 approaches zero, but that motion could still be driven by bathymetric variations. Assuming a separable wave solution of the form: $ = Z(z)ei(kX+ly"at) (41) leads to the following Sturm-Liouville problem for the verti- cal dependence: '^ Z I + [k2 + I2 + & 6P z/z L Z = 0 (42) 30 Substituting from eq . (7) for p and 6 , using the fact that the Brunt-Vaisala frequency, N, is given by converting back to a dimensional form results in: , 2 g Hs pq 8z and o s u N' k2 + 2 + Bk a Z = 0 (43) B. EQUATION USING ARCTIC DFNSITY PROFILE 1 . Basic State Description In this initial formulation zero mean flow was assumed, thus the initial density distribution was assumed to be a function of depth only. The Arctic under-ice density profile was approximated by a homogeneous upper layer from the surface to z = -h, below which density increased follow- ing an exponential curve. In this lower layer density increases such that it assymptotically approaches a constant, p , , at deeper depths. This distribution is similar to that proposed by Reid [1948] , in his equatorial model, and is given analytically by: Ks "n ps = pd " (pd"ph)e (0 > z > -h) (44) l/c(z+h) (z < -h) which is illustrated in Figure 8 for various values of c and h (referred to hereafter as layer depth (h) and e-fold depth (c) ) . Actual density values taken during the collection of velocity data are also shown on Figure 8 for comparison. As can be seen, the best data fit is with an e-fold depth of between 100 and 120 meters and a layer depth of 30 meters. 31 100 200 — a a 500 600 300- 400 — X - Gait [1967] 700 Figure 8. Depth Distribution of Analytic Density Function Shov.'ing Relation to Data. 32 2 . Use of Actual Density Distribution in Equations The density distribution given by eq. (44) was sub- stituted into eq. (43) to produce a vertical normal mode equation which would apply to this specific density structure. To show this, eq . (43) is first expanded in terms of p and its derivatives : o -^— Z + 2Z - p zz z L. S p p s*s zz Ups ) z ZJ k2.+l2 + ^ a Z = 0 (45) (46) The derivatives of p are found from eq . (4 4) as "^d"^ l/c(z+h) p_ = e (pd ph} l/c(z+h) = 2 e zz c which can then be substituted into eq. (44) to yield, after some rearrangement of terms : Z zz "2(pd-ph} l/c(z+h) M 1 e + — p C CJ f c o 2~ (pd-ph' k2 + I2 + to el/c(z+h)z (47) This can be simplified to: zz Ael/c(z+h) + z + XeVc(z+h) . Q (48) where A = 2(pfl-Ph) P c s (49) 33 B = (49) 2 f cp o Ks (pd-ph} 2 2 k + 1 + Pk The solutions to eq. (48) will then be a set of eigen- functions describing the vertical distribution of the pressure function, $ . It is noted that this equation only applies to the region below the layer depth. In the upper homogeneous region horizontal velocities will be constant since there is no stratification. Vertical boundary conditions on the system are that w, the vertical velocity, is equal to zero at both the surface (ice) and bottom boundaries. In addition, velocities are assumed to be continuous across the internal boundary at z = -h. If these boundary conditions are met and a solution to eq. (48) exists there will be solutions for different discrete values of A, the eigenvalues for this equation. A value of either k or a may then be assumed and used with the eigenvalues to solve for the corresponding values of a or k for each mode. In order to actually solve for both k and a some initial distribution for can be written as: (J) = Z (z) cos (kx-at) From eq. (36) w can be written as: yps w = r Z sin (kx-at) 0 z p which becomes : f w = — y Z sin (kx-at) (52) IT z when the terms are dimensionalized and proper substitutions made for y and 6 . Assuming a value for k is known, an expression for a can then be derived from eq . (49) : 35 0 = —2 (53) f cp X „ g(Pd-Ph) This is similar to the equation for the frequency of a baro- tropic Rossby wave, a, = $/k , v;ith the additional term in the denominator being due to the stratification. 36 III. RESULTS OF WAVE SOLUTION AFP COMPARISON WITH ACTUAL DATA A. NUMERICAL SOLUTION OF THE VERTICAL MODE EQUATION The vertical mode equation (48) was rewritten using a normalized depth variable, 0 >_ z ' >_ -1.0, as: I - D(AeDB(z'+h/D)+B)Z- -X'eDB(z'+h/D»Z=0 (54) zz z 2 where X1 = D X and D = depth of water. Equation (54) was solved numerically on the Naval Post- graduate School IBM/360 computer. The solution was found using a Milne predictor-corrector method, according to an algorithm described by Nielsen [1964]. The depth range was divided into 2000 increments and then the Milne method used from the bottom, z1 = -1.0, upward to the normalized layer depth, z1 = - h/D. The Milne method section of the numerical solution was then tested by using it to solve a zero order Bessel equation which had a known analytic solution. Results of this test were satisfactorily checked with a table of Bessel functions. In the relatively thin homogeneous region above the layer depth, the slope of the vertical velocity was. held at its layer depth value, as was the value of the vertical structure function, Z. The output of the numerical scheme was a set of 2000 values of Z and Z from z1 = -1.0 to z' = 0, for a given z X*. The vertical distribution of horizontal velocity is proportional to Z as was shown by eq . (50). From eq . (52) it 37 can be seen that the vertical dependence of vertical velocity, 2 w, is proportional to 7, /N . The numerical scheme could thus be used to give the vertical dependence of both horizontal and vertical velocities. The bottom boundary condition, w = 0 at z' = -1.0, was met by setting Z (-1.0) = 0. This was valid since w is 2 2 proportional to Z /N and since N (-1.0), although quite small, is not identically zero. With the bottom boundary condition thus fixed, the condition w = 0 at z ' = 0 could only be met for certain values of X ' . These values are the eigenvalues corresponding to the normal modes of oscillation of the ocean with the given vertical density distribution. In order to find the correct eigenvalues an initial X ' was assumed. An iterative process was then begun in which the numerical method described earlier was used to arrive at a value of w(0) for each new value of X ' . The initial value chosen for X' was then increased in steps until the value of. w(0) changed sign. At that point a Regula-Falsi method (eq. (55) below) was used to find the next X'. X' w (0) - X'w , (0) \ > - n~1 n n n-1 , r c- x n+1 w (0) - w , (0) KDD> n n-1 This was continued until the value of w(0) was sufficiently close to zero. The value of X ' at that point was then con- sidered an eigenvalue. The process was continued by increas- ing the value of X ' in steps until the next eigenvalue was found in the same manner. 38 B. RESULTS OF THEORETICAL DFVELOPMFNT Figures 9 and 10 show the vertical distribution of horizontal and vertical velocities which result from the numerical solution of eq. (54) for a given density distribu- tion (h=30, c=120) . The first four normal modes are shown in each figure. Velocity values plotted are not absolute values but are relative to the horizontal velocity at the bottom for each mode. The ratio of w(z) to V(D) was computed using a and k values corresponding to a wavelength of 30 km. As can be seen on these graphs, the vertical velocity would be insignificant for a horizontal velocity on the order of cm/sec or larger. This means the only measurable motion would be in the horizontal direction. Figure 11 shows the effect of a change in the density distribution layer depth on the horizontal velocity distribu- tion. The second mode only is shown since effects on other modes are similar. The values chosen, 20 and 50 meters, represent the normal seasonal variation in the thickness of the under-ice layer. As can be seen, very little change in the velocity distribution results. Figure 12 shows the effect on the second mode of a change in the curvature (e-fold depth, c) of the below layer density distribution. This shows that as the e-fold depth increases (decrease in curvature) , the velocity peak is less marked and occurs at a deeper depth. Frequencies were solved for by using eq . (53) with an assumed horizontal wavelength of 30 km. Results were on the 39 Horizonlal Velocify [ vfz) /v ^30003) 1000-t 1200 hs30m co120m Figure 9. Vertical Structure of Theoretical Horizontal Velocity Patterns for Modes; 1 through 4. 4 0 -1*10 r»4 Vertical Velocity [w(0/ v C3000:j) x-4 3000 Figure 10. Vertical Structure of Theoretical Vertical Velocity Patterns for Modes 1 through 4. 41 Figure 11. Effect of Layer Depth Change on Vertical Structure 4 2 Horizontal Velocity £ v Cz)/vC3000) ] -15 -10 -5 0 5 Figure 12. Effect of Varying e-fold Depth in Density Distribution on Vertical Structure. — 8 — 1 order of 10 sec or less, which corresponds to periods of tens of years. These values indicate the current patterns are essentially stationary. As can he seen from eq . (53) either a decrease in wavelength (increase in k) or an increase in mode (higher A) will cause the frequency to decrease. Thus the limiting, highest possible frequency for a given density distribution and horizontal wavenumber is for A = 0, a barotropic Rossby wave mode, and is given by: b k For a wavelength of 30 km this corresponds to a period of about 10 years. The horizontal distribution of velocities was arrived at by assuming a length of 15 km to the first zero of the sin(kx-ot) function of eq. (50) and forming the product of this trigonometric function and the Z (z) function at a con- stant depth. This is shown in Figure 13. C. COMPARISON WITH OBSERVED DATA A comparison of the theoretical normal modes plotted in Figure 9 with the observed data in Figure 2 shows that the second normal mode most nearly matches the measured vertical distribution of horizontal velocity. The only real discre- pancy between the two is the relative maximum in the theo- retical distribution at the surface. This appears to be an inherent characteristic of the theoretical solution and is seen in each of the normal modes. One possible explanation is that these theoretical maxima were caused by treatment of 44 % I I 1 t. L/2 Figure 13. Horizontal Distribution of Theoretical Horizontal Velocities. 45 the under-ice surface as flat and frictionless in the develop- ment. The ice actually has a rough surface capable of dis- sipating some energy of relatively slow moving currents. A second explanation is that more than one mode exists with a progressively decreasing amount of energy in succes- sively higher modes. Thus if v (z) represents the horizontal velocity distribution for a given mode, then: oc v (z) = S a v (z ) , , n n n=l where a > a ,,, would describe the actual velocity distri- n n+1 J bution. Assuming further, as a simplistic first approach, that the only significant a are for n=l and n=2 (i.e., only the first two modes contain significant energy) the linear combination of modes 1 and 2 was plotted in Figure 14. The a, and a„ were sized to give zero velocity at the surface, meaning a~ was less than a, . Figure 14 can be seen to be closer to the observed data than the second mode alone. Obviously a more rigorous solution for the coefficients could be used to get an even closer fit to the data. This, however would not be physically meaningful. The purpose of this initial investigation was merely to show that the theoretical development would allow motion similar to what has been observed. This has been shown. The next phase of the inves- tigation, which will be discussed in Section IV of this thesis, is to study the initial value problem of how this specific ocean system responds to the initial perturbations in its density structure caused by a newly forming lead (i.e., which normal modes, if any, are excited). 46 Horizontal Velocity [ vCz)/v(3000)3 0 10 20 Figure 14. Linear Sum of Fodes 1 and ? for h = 30, c - 120. In order to compare real and theoretical distributions more closely, a vertical cross-section through the origin was drawn showing theoretical horizontal velocity values, Figure 15. This was done for the second normal mode for a density distribution having h=30 and c=120, and for the simple horizontal distribution of Figure 13. Thus Figure 15 could be compared with Figure 4, which is drawn for real current data having a similar density distribution. A similarity can be seen in the general flow direction and in the shape of the isolines. Apart from the relative smoothness of the theoretical plot, the only significant difference is in the upper homogeneous layer, which has been discussed previously. Figure 16 is the same cross-section drawn for the summation of the first two normal modes. As can be seen, the cell structure is closer to the observed data of Figure 6 than was the second mode alone. 48 800 Horizontal Distance L/2 ~sr * ^sinCku>! vC3000) c = 120m I, -. 30 m Figure 15. Vertical Cross-section of Horizontal Velocities, Mode 2. 40 200 a. 400 v a 600 800- Horizonfal Distance L/2 - , —J c = 120m h« 30m Figure 16. Vertical Cross-section of Horizontal Velocities, Sum of Modes 1 and 2. 5 0 IV. DEVELOPMENT OF THE INITIAL VALUE PROBLEM A. GENERAL METHOD The problem addressed here is whether the density field perturbations which occur at a newly forming Arctic lead can force the theoretically developed oscillations to occur and, if so, whether the forced oscillations are consistent with observed data. The solution to this problem lies in formulating expli- citly the initial effect of the lead-induced perturbations on the pressure field. This effect can then be used as an initial condition on the pressure function, cf> , in order to fully solve the Sturm-Liouville problem developed in Sections II and III. It was shown that the general solution to that problem was a set of eigenfunctions which gave the initial response of the Arctic Ocean to a baroclinic perturbation. The linear combination of these functions can be written as: oc (x,z,t) = 2 ' S A B Z (z)cos(k x-a t) (56) Y ., n m n m ran n=l m=l so that at time zero: oo o: (x,z) = 2 £ A B Z (z)cos(k x) (57) Yo , , n m n m n=l m=l It must be noted that the summation of solutions indicated by eq . (56) is not necessarily a solution itself. It would only constitute a solution if the initial equations were actually 51 linear, which they are not. It appears reasonable to assume however that the non-linear contribution is small compared to the linear so that this summation approximates an actual solution. A method will now be shown whereby the values of A ,B ,k and a can be found based on the initial specifica- n m m mn L tion of the ^/^ Figure 17. Perturbation Density Due to Salt Rejection, at Tine Zero. " ■■-■" x Figure 18. Vertical Cross-section of Perturbation Density at Time Zero. 53 -;.- z initial overpressure Figure 19. Perturbation Pressure at Time Zero. 54 X *^ X ♦0(x,z) = - ~ I ' P (58) 0 < z < z - - p X I < X = 1 I - p z > z p; =0 ( Ixl > x ) 1 ' p' C. SOLUTION OF THE PROBLEM The initial condition, cf> (x,z), is separable and can be written as x ) (59) =1 ( Ixl < x ) g(z) = -z/h (z > -h) (60) =1 (z < -h) Thus, coefficients can be solved for separately from: f (x) = Z B cos (k x) , m m m=l (61) oo g ( z ) = E A Z ( z ) t n n n=l Considering B first, these are the coefficients for a Fourier cosine series, which on the interval (0,L) are given by: x -i P i L x B^ = - / f (x)dx + =- f f (x)dx = -£■ (62) ° L 0 L x L x p Bm - | / P f (x)cos(~^)dx (111=1,2,3...) = ~ sin(kmxp) (m=l,2,3...) (63) 55 where k = miT/L is a distance sufficiently far removed from m J the lead so that 4> can be assumed to have a zero value. The values of B correspond to the zeroth mode, which is o ^ a barotropic mode that may exist but is not of interest in the present discussion. Figure 20 is a plot of the approx- imation to f (x) formed by evaluating eg. (61) for m = 1 through m = 40. It can be seen to approximate f (x) fairly closely. Figure 21 is a plot of the horizontal distribution of v for these same modes. For each k thus found, and each eigenvalue corresponding to the Z (z) vertical structure functions, there can be found n a value of a mn The values of A are given by: n ^ J A = n -D / g(z)Z (z)dz 0 -D (n=l , 2 , 3 . . . ) (64) 2 /. [Z (z)]"dz 0 n for a general series expansion, assuming the Z (z) are orthogonal functions with respect to a weighting function of one. The denominator can be normalized for all modes by finding a factor C such that: n "D 2 / [C Z (z)rdz = 1 o n n or : n "D 2 ' / (?, u)rdz - o (65) Thus A becomes: n 56 too 10 T -5 xp=1 L=20 i«nn 10 -'-•x Figure 20. Fourier Series Fxpansion of Initial Horizontal Pressure Distribution. vCx) "> -x Figure 21. Fourier Series Fxpansion of Horizontal Velocity Distribution in ,x-direction . 57 -D A = / g(z)C Z (z)dz (n=l,2,3...) (66) n -, n n for the normalized distributions. The computation of A can be further simplified by writing g(z) as the sum of an average and a variable part. Since D » h, the average part is essentially: g (z) - 1. ^avg The only variable part is then: g.ra„(z) = 1 + z/h (z > -h) (67) Thus, eq. (64) can be written: -D -h A = f C Z (z)dz + / (l+z/h)C Z (z)dz (68) n ' nn ' ' ' n n But, the first integral in eq . (68) is equal to zero based on the original continuity relationship used in deriving Z(z). From Figure 9 it can be seen that each Z (z) is constant for a n depths above -h. Thus eq . (68) finally becomes: A = - h/2 [a C ] (n=l,2,3...) (69) n n n These values were solved for using the IBM/360 computer and were then used to form the series summation given by eq. (61) for the first six modes. The results are shown in Figure 22. The approximation is not as close as the horizontal case due to using only six modes. This was a limit imposed by the numerical method. A complete set of values for A , B , k , and a is then ^ n m m mn available. Values for the distribution having c=120, h=30, -1 500 a CI 1000 sc 9var^ A„Z„C«> Figure 22. Series Expansion of Vertical Structure Function. 59 are given in Tables I, II, and III. Values of L=20 km and r =1 km were assumed. P The solution for 4> given by eq . (56) will then begin with a distribution as given by eq . (57), which is the ini- tial pressure distribution. Due to the extremely small values of a there will be no significant movement with time of the overall wave system. Phase velocities (a/k) are on _2 the order of 10 cm/sec. The distribution of horizontal velocities however, is proportional to: oo oo v(x,z,t) = E 7, A B k Z (z)sin(k x-c t) (70) ' ' ' , ,nmmn m mn n=l m=l The initial horizontal and vertical distributions of v(x,z,t) were shown in Figure 21 and Figure 22. TABLF I VALUFS OF o (sec~1xl0~9) AND T (yrs) FOR mn mn J. SELECTED VERTICAL AND HORIZONTAL MODES m 10 1 9.8 3 5.84 3.11 0 mn 20.0 33.9 63.4 T mn n 3 2.58 4.38 2.85 a mn 76.9 45.2 69.3 T mn 5 1.11 3.0 2 2.49 a mn 178. 65.6 77.0 T mn 60 TABLE II VALUES OF k AND B FOR HORIZONTAL m m MODES 1 THROUGH 10 m _ 1 _ f- km(cm xlO ) Bm 1 1.56 .0995 2 3.14 .0983 3 4.70 • .0963 4 6.28 .0935 5 , 7.85 .0900 6 9.40 .0858 7 10.1 .0810 8 12.5 .0756 9 14.1 .0699 10 15.6 .0637 TABLE III VALUES OF X AND A FOR VERTICAL n n MODES 1 THROUGH 6 n n n 1 +.0237 9.18xl0-5 2 -.0218 4.66xl0~4 3 +.0193 1.14xl0~3 4 -.0169 2.14xl0~3 5 +.0150 3.47xl0~3 6 -.0134 5.15xl0~3 61 V. CONCLUSIONS The stated objective of this study was to model certain observed Arctic Ocean internal current patterns through an analytic study of the baroclinic response of the ocean to lead formation and refreezing. The problem was then divided into two major phases: (1) finding the baroclinic internal normal modes of the Arctic Ocean; and (2) investigation of the forcing of these modes by the effects of lead formation. The normal modes of response found as a result of the theory developed in Section II showed a velocity pattern which was essentially stationary with no significant vertical velocity. However, horizontal velocities of the magnitude measured in the Arctic could exist within these stationary patterns. It was noted that the steepness of the pycnocline (measured by the parameter, c) had a direct and significant effect on the depth of maximum velocity, whereas the effect of a variation in layer depth was relatively minor. The apparent consistency of the cell structure seen in Figure 16, the theoretical vertical cross-section of horizontal veloci- ties , with the observed data lends strength to the plaus- ibility of the developed theory. The fact that a simple linear combination of the first two theoretical normal modes, Figure 14, qualitatively matches measured velocity profiles is of particular note. It shows the physics used to develop the theory to be 62 consistent with actual data, in that observed current patterns appear to be part of the set of allowed solutions. It has not yet been shown conclusively that the perturbations at a lead will force these modes. The study done here is a first step which looked at the dynamics of a single Arctic lead. The .theory developed has been shown to be dynamically possible but the study has not yet looked at the energy involved in the system. It is possible that an interaction between leads or between lead(s) and the mean flow is needed to provide the energy required for the observed current magnitudes. The very small values of frequency, corresponding to periods of greater than 20 years, are attributable to the magnitude of 3 in the Arctic. The vorticity equation derived, eq. (38), is a non-linear time dependent parabolic equation. However, the results found in the initial value problem studied in Section IV of this thesis indicate that the time dependence of the equation is not significant since the solu- tions were essentially stationary. Under the same geometry the problem could be reformulated as an elliptic , linear equation forced with time dependent boundary conditions. An advantage of this type formulation would be the ability to write a separable differential equation in cylindrical coordinates , which was not possible using the development of Section II. Motion approximating the vortex-like patterns observed by Coachman and Newton [197 3] could then be more easily described. 63 The plane wave solution for velocities developed in this thesis is a result of the formulation of the problem in rectangular coordinates with an infinite length lead. This worked well in developing the theory and resulted in a real- istic vertical section of velocity structure. However, actual leads are of finite length which could be expected to induce curvature in the velocity patterns. 64 VI . RECOMMENDATIONS FOR FURTHFR RFSFARCII The consistency found with observed data warrants the continuation of this study using the theory developed in this thesis as a basis. The following recommendations for further research are made: 1. the formulation of the initial value "study as an elliptic problem, as discussed in section V; 2. the use of vertical stretching caused by bathymetry, or forced convection at a lead, rather than the 3 effect as the driving mechanism in the developed vorticity equation; 3. the investigation of the question of whether these velocity patterns have any influence on the stability of the larger scale mean flow, including the effects of baroclinic instabilities ; 4. a study of the energetics of the system. 65 BIBLIOGRAPHY Coachman, L. K. , "Production of Supercooled Water Durina Sea Ice Formation," Proc . Symp . Arctic Heat Budget and Atmos- pheric Circulation , The Rand Corporation, RM-5233-NSF, 497-529, 1966. Coachman, L. K. and J. L. Newton, Personal Communication, 1973. Foster, T . D. , "Haline Convection Induced by the Freezing of Sea Water," Jour. Geophys . Res., v. 73, 193 3- 1938, March 1968. Foster, T. D., "Haline Convection in Polynyas and Leads," Jour. Phys. Ocean. , v. 2, 462-469, October 1972. Gait, J. A., Current Measurements in the Canadian Basin of the Arctic Ocean, University of Washington Department of Oceanography Technical Report 184, 1967. Gait, J. A. , "A Numerical Investigation of Arctic Ocean Dynamics," Jour. Phys. Ocean., In print, October 1973a. Gait, J. A., Unpublished Manuscript, 1973b. Hildebrand, F. B. , Advanced Calculus for Applications, Prentice-Hall, Inc., 1962. Kuo, H. L., "Dynamics of Quasigeostrophic Flows and Insta- bility Theory," Advances in Applied Mechanics, v. 13, 248-330, Academic Press, 1973. Nielsen, Kaj L., Methods in Numerical Analysis, MacMillan, 1964. Pedlosky, J., "The Stability of Currents in the Atmosphere and Ocean: Part I," Jour. Atmos . Sci. , v. 21, 201-219, March 1964. Pedlosky, J., "Geophysical Fluid Dynamics," Mathematical Problems in the Geophysical Sciences, 49-60, American Mathematical Society, 1971. Piascek, S. A., Numerical Studies of Convection Currents at the Ocean's Surface Caused by Fvaporation , Radiation and Atmospheric Cooling , Final Technical Report for ONR Contract N0"0014-67-A-0242-003 , 19 70. 66 Reid, R. O., "A Model of the Vertical Structure of Mass in Equatorial Wind-driven Currents of a Baroclinic Ocean," Jour. Marine Res., v. 7, 304-312, November 1948. Schaus , R. H. , A Thermodynamic Model of a Central Arctic Open Lead, M.S. Thesis, Naval Postgraduate School, 1971. Semtner, A. J., Personal Communication, 1973. 67 INITIAL DISTRIBUTION LIST No. Copies 1. Defense Documentation Center Cameron Station Alexandria, Virginia 22314 2. Oceanographer of the Navy Hoffman Building No. 2 2461 Eisenhower Avenue Alexandria, Virginia 22314 3. Office of Naval Research Code 4 80 D Arlington, Virginia 22217 4 . Naval Postgraduate School Code 58 Monterey, California 93940 5. Library, Code 0 212 Naval Postgraduate School Monterey, California 93940 6. Dr. R. E. Stevenson Scientific Liaison Office Scripps Institution of Oceanography La Jolla, California 92037 7. Dr. Jerry A. Gait NOAA Pacific Marine Env. Lab. University of Washington WB-10 Seattle, Washington 98105 8. LCDR Richard I. Itkin, USN 96 Bonaire Circle Suffern, New York 10901 9. R. K. McGregor Director Arctic Program, Code 415 Office of Naval Research Arlington, Virginia 22217 10. AIDJEX Division of Marine Resources University of Washington Seattle, Washington 9 8105 68 11. Assoc. Professor R. T. Williams Department of Meteorology Naval Postgraduate School Monterey, California 93940 12. Professor L. K. Coachman Department of Oceanography University of Washington, WB-10 Seattle, Washington 98105 6^ UNCLASSIFIED SECURITY CLASSIFICATION OF THIS PAGE (Whan Dale Fntered) REPORT DOCUMENTATION PAGE READ INSTRUCTIONS BEFORE COMPLETING FORM 1. REPORT NUMBER 2. GOVT ACCESSION NO. 3. RECIPIENT'S CATALOG NUMBER 4. TITLE (and Subtitle) AN ANALYTIC MODEL OF OBSERVED SMALL SCALE BAROCLINIC CURRENTS IN TEE ARCTIC OCEAN 5. TYPE OF REPORT & PERIOD COVERED Master's Thesis; September 197 3 6. PERFORMING ORG. REPORT NUMBER 7. AUTHORfs; Richard Ivan Itkin; LCDR, USN 0. CONTRACT OR GRANT NUMBER(»J 9. PERFORMING ORGANIZATION NAME AND ADDRESS Naval Postgraduate School Monterey, California 93940 10. PROGRAM ELEMENT. PROJECT, TASK AREA & WORK UNIT NUMBERS II. CONTROLLING OFFICE NAME AND ADDRESS Naval Postgraduate School Monterey, California 93940 12. REPORT DATE September 19 7 3 13. NUMBER OF PAGES 71 H. MONITORING AGENCY NAME A ADDRESSf// different from ConttolUn6 Office) IS. SECURITY CLASS, (of this report) Unclassified 15». DECLASSIFI CATION/ DOWN GRADING SCHEDULE 16. DISTRIBUTION ST ATEMEN T (of thte Report) Approved for public release; distribution unlimited. 17. DISTRIBUTION STATEMENT (of tho abstract entered In Block 20, If different from Rnport) 18. SUPPLEMENTARY NOTES 19. KEY WORDS (Continue on reverse eldo If necessury and Identify by block nu»n Dr.te Entered) 71 '^H'15 5 4593C ,ode\ of e in ^V* ltk»n &n < obse , . r cui • — *ed -S*acurr< tbe Thesis 145930 175 Itkin c.l An analytic model of observed small scale baroclinic currents in the Arctic Ocean. thesl75 An analytic model of nh mZmm^ fiscal 3 2768 002 101 So