1 A NORMAL MODE MODEL FOR ESTIMATING LOW-FREQUENCY ACOUSTIC TRANSMISSION LOSS IN THE DEEP OCEAN Ki rk Eden Evans OUOLEY K*l MAVD.X. POS NAVAL POSTGRADUATE SCHOOL Monterey, California THESIS A 1 FORMAL MODE MODEL FOR ESTIMATING LOW -FREQUENCY ACOUSTIC TRANSMISSION LOSS IN THE DEEP OCEAN by Kirk Eden Evans September 1975 Th esis Advisor: Alan B. Coppens Approved for public release; distribution unlimited. U173522 SECURITY CLASSIFICATION OF THIS PAGE (Whan Data Enl.r.d, REPORT DOCUMENTATION PAGE READ INSTRUCTIONS BEFORE COMPLETING FORM 1. REPORT NUMBER 2. GOVT ACCESSION NO 3. RECIPIENT'S CATALOG NUMBER 4. TITLE (and Subtitle) 3. TYPE OF REPORT & PERIOD COVERED A Normal Mode Model for Estimating Low-Frequency Acoustic Transmission Loss in the Deep Ocean Ph.D. Thesis September ]975 6. PERFORMING ORG. REPORT NUMBER 7. AUTHORS Kirk Eden Evans 8. CONTRACT OR GRANT NUMBERfaj 9. PERFORMING ORGANIZATION NAME ANO ADDRESS Naval Postgraduate School Monterey, California 93940 10. PROGRAM ELEMENT, PROJECT, TASK AREA a WORK UNIT NUMBERS II. CONTROLLING OFFICE NAME ANO ADDRESS Naval Postgraduate School Monterey, California 93940 12. REPORT DATE September 1975 13. NUMBER OF PAGES 178 14. MONITORING AGENCY NAME ft AOORESSf// dlllarant from Controlling Ofllca) Naval Postgraduate School Monterey, California 93940 IS. SECURITY CLASS, (ol thta tiport) Unclassified I5«. DECLASSIFI CATION/ DOWN GRADING SCHEDULE 16. DISTRIBUTION STATEMENT (ol thla Raport) Approved for Public Release; Distribution Unlimited. 17. DISTRIBUTION STATEMENT (ol tha abatract antarad In Block 30, II dlltarant from Raport) 18. SUPPLEMENTARY NOTES 19. KEY WORDS (Contlnua on ravaraa alda II nacaaaary and Idantlty by block nuotbar) Underwater sound Normal Modes WKB SOFAR PARKA Transmission Loss Underwater Acoustics Low- frequency LOFAR Ray-Mode Sound Channel Acoustic Propagation 20. ABSTRACT (Contlnua on ravaraa alda it nacaaaary and Idantlty by block numbar) This work describes a computer model (QMODE) which uses the normal mode method for the estimation of low-frequency long-range acoustic transmission loss in the deep ocean. For computational speed, the model uses the Wentzel-Kramer- Brillouin (WKB) solutions for the phase speeds (eigenvalues) and modes (eigenf unctions) . The WKB solutions are extended DD , It?-,, 1473 (Page 1) EDITION OF 1 NOV 69 IS OBSOLETE S/N 0 102-014- 660 1 | SECURITY CLASSIFICATION OF THIS PAGe (Whan Data gntarad) SbCUKITY CLASSIFICATION OF THIS PtGEfW^an Dvtm En f rod 20. ABSTRACT (cont) to consider the effects of the surface and bottom boundaries. The eigenvalues are initially estimated by a least-squares polynomial fit through sample results of the WKB characteristic equation. The effects of range dependence in the sound speed profile are simulated through the use of the adiabatic assumption, which has been extended to include the case of profiles with multiple sound channels. Results of the model agree generally with the transmission loss observed in two long-range acoustic experiments. DD Form 1473 . 1 Jan 73 S/N 0102-014-6601 SECURITY CLASSIFICATION OF THIS PAOEOWiw Data Enfrtd) A NORMAL MODE MODEL FOR ESTIMATING LOW-FREQUENCY ACOUSTIC TRANSMISSION LOSS IN THE DEEP OCEAN by Kirk Eden Evans Lieutenant Commander, united States Navy 5. A., Miami University. 1966 M.S., Naval Postgraduate School, 1973 Submitted in partial fulfillment of the requirements for the degreeof DOCTOR OF PHILOSOPHY from the NAVAL POSTGRADUATE SCHOOL September 1975 ■ ABSTRACT This work describes a computer model (QMODE) which uses the normal mode method for the estimation of low-f reguency long-range acoustic transmission loss in the deep ocean. Hentzel-Kramer-Brillouin (HKB) solutions for the phase speeds (eigenvalues) and modes (eigenf unctions) . The WKB solutions are extended to consider the effects of the surface and bottom boundaries. The eigenvalues are initially estimated by a least-squares polynomial fit through sample results of the IfKB characteristic- eguation. The effects of range dependence in the sound speed profile are simulated through the use of the adiabatic assumption, which has been extended to include the case of profiles with multiple sound channels. Results of the model agree generally with the transmission loss observed in two long-range acoustic experiments. TABLE OF CONTENTS LIST 0? TABLES 9 LIST OF FIGURES 10 I. INTRODUCTION.- , 13 A. DEVELOPMENT OF TEE NORMAL MODE TECHNIQUE 14 1. Early Theory 14 2. Numerical Methods 19 3. Current Models 21 B. OPERATIONAL MODELS 22 C. SCOPE OF THIS WORK 23 1. Objectives 23 2. Departure from Previous Works- 23 D. REFERENCES 25 II. THEORY: HOMOGENEOUS WAVE EQUATION 28 A. HOMOGENEOUS WAVE EQUATION 31 1. The Helmholtz Equation -.. 32 2. The Hankel Transform 32 3. Nature of the Solution - -.. 33 B. THE WKB METHOD 35 1. The WKB Solution 35 2. WKB Validity Requirements 37 3. Properties of the WKB Solutions 38 4. WKB Characteristic Equation.., 38 C. REFERENCES „ 41 III. THEORY: AIRY SOLUTION AND CONNECTION FORMULAE 42 A. THE AIRY SOLUTION , 42 1. Airy Solution for Small Argument 43 2. Airy Validity Factor - 43 3. Airy Solution for Large Argument 44 B. CONNECTION FORMULAE 46 1. WKB Oscillatory and Airy Connection...... 46 2. WKB Exponential and Airy Connection 47 3. h'KB Connection Formulae . 47 C. BARRIER CONNECTION FORMULAE 47 D. REFERENCES 51 IV. BOUNDARY CONDITIONS 52 A. THE SURFACE BOUNDARY 52 1. Case I .„ 52 2- Case II 54 3. Case III . 57 4. General Cases. . 57 B. BOTTOM BOUNDARY CONDITION 62 1. Case 1 64 2. Case II 66 3. Case III „ ., 68 4. The General Cases .. 70 C. BOUNDARY CONDITION TEST 70 D. REFERENCES «... 73 V. THEORY: THE INHOMOSENEOUS EQUATION 80 A. SOLUTION OF THE INHOMOGENSOUS EQUATION , . 80 1. The Green's Function f.... 80 2. Evaluation of the Integral....- 82 3. Normalization 36 B. THE Q FUNCTION 87 C. REFERENCES 87 VI. THE MODEL .• ... 88 A. MODE FAMILIES 89 1. Definition of Families 89 2. Family Types, 92 3. Relationship of Q between Families 95 a. Single Channel 95 b. Double Channel 95 B. EIGENVALUE ESTIMATION 99 1. The Eigenvalue Search 99 2. Selection of the Estimator - 100 a. Cubic Spline..... 100 b. Least Squares Polynomial. 101 C. CALCULATION OP THE EIGENFUNCTION 102 1. 'Assumed Solutions 102 2. Spliced Solution 103 3. Eigenf unction Normalization..- 105 D. TRANSMISSION LOSS 105 1. Coherent Transmission Loss „... 106 2. Averaged Transmission Loss.... 107 3. Incoherent Transmission Loss.- 108 E. MULTIPLE PROFILE PROPAGATION., 108 1. Adiabatic Assumption .„ 109 2. Adiabatic Transition Range.... 109 3. Transmission Loss 110 4. Mode Matching 111 a. Same Channel 114 b. Cross-Channel Transitions 114 F. BOTTOM MODELS 116 1. Full Fluid Bottom 116 2. Rigid Bottom 117 3. Modified Fluid and Rigid Bottoms 118 4. Bottom and Surface Roughness... 118 G. REFERENCES 119 VII. THE PROGRAM „ 120 A. MAIN PROGRAM .. 120 B. SUBPROGRAMS 123 1. BOTSET 123 2. BCSET 124 3. FINFAM 124 4. ORTPOL.. 125 5. EIGFUN 125 6. SPLICE 126 7. PLOSS 127 8. MATCH 128 9. INTEGR 129 10. QFIND 129 C. REFERENCES 129 VIII. RESULTS AND CONCLUSIONS 130 A. ANALYTIC TESTS 130 1. Isovelocity Channel 130 2. Constant Gradient Channel.... 131 B. AESD WORKSHOP TESTS 131 1. Comparison of Moles 132 2. Transmission Loss 132 C. COMPARISONS WITH DATA.. 133 1. PARKA IIA Data 133 a. Run 1.... • 133 b. Run 2....o 134 c. Run 3....- 134 d. Run 4 134 2. Antigua to Grand Banks - 135 D. EFFECT OF THE BOTTOM 136 E. CONCLUSIONS 137 F. REFERENCES 139 INITIAL DISTRIBUTION LIST ..., 175 LIST OF TABLES A. AESD 1 Profile and Run Times 77 B. AESD 2 Profile and Run Times 78 C AESD 3 Profile and Run Times 79 LIST OF FIGURES 1. Path of Bessel Integrals in the Complex Plane 16 2. Sound Speed Profile 29 3. The BKB Solutions 40 4. Validity Criteria as a Function of Scaled Depth..... 45 5. Double Sound Channels and Barrier - 49 6. Surface Boundary Condition - Case I 53 7. Surface Boundary Condition - Case II 55 8. Surface Boundary Condition - Case III... 58 9. Surface Boundary Condition - -9 vs x.... 59 u 10. Typical Sound Speed Profile near the Surface 61 11. Surface Boundary Condition - Q vs c 63 u p 12. Bottom Boundary Condition - Case 1 65 13. Bottom Boundary Condition - Case II..... 67 14. Bottom Boundary Condition - Case III.... 69 15. Bottom Boundary Condition - Q vs c 71 1 P 16. ASSD 1, e and % vs c 74 1 up 17. AESD 2, e and 9 vs c 75 1 up 10 18. AESD 3, 6 and Q vs c , 76 1 a p 19. Path of the Hankel Integrals in the Complex k Plane - „ 84 r 20. Bilinear Profile Families 91 21. Family Classifications I . c. 93 22. Family Classifications II 94 23. Double Channel with Barrier, 96 24. Ill-Conditioned Cubic Spline 101 25. The Splicing Technique ....104 26. Two Near Symmetric Channels with Modes.- 112 27. Cross Channel Transitions...., 113 28. Eigenvalue Test: Isovelocity Channel...- 140 29. Constant Gradient Channel -- 141 30. Eigenvalue Test: Constant Gradient Channel 142 31. AESD 1 - Sound Speed Profile 143 32. AESD 2 - Sound Speed Profile 144 33. AESD 3 - Sound Speed Profile 145 34. AESD 3 Modes 146 35. AESD 2 Modes 11 through 15 147 36. AESD 2 Modes 16 through 20 148 37. AESD 1 - QMODE Results. 149 38. AESD 1 - FACT Results ....150 39. AESD 1 - NOL Results 151 40. AESD 1 - Parabolic Equation Results...... 152 11 41. AESD 2 - QMODE Results.- .. 153 42. AESD 2 - FACT Results . 154 43. AESD 2 - NOL Results.-. . ....155 44. AESD 2 - Parabolic Eguation Results 156 45. AESD 3 - QMODE Results.. , 157 46. AESD 3 - NOL Results 158 47. AESD 3 - Parabolic Equation Results 159 48. AESD 1 - 2.5 nm Averaging 160 49. AESD 1 - 5.0 nm Averaging f 161 50. AESD 1 - 7.5 nm Averaging. . « - -...162 51. AESD 1 - 10 nm Averaging 163 52. PARKA IIA Run 1 Results 164 53. PARKA IIA Run 2 Outbound 165 54. PARKA IIA Run 2 Inbound 166 55. PARKA IIA Run 3 Results '...167 56. PARKA IIA Run 4 Results ....168 57. Antigua to Grand Banks - Bilinear Profile 169 58. 13.89 Hz Bilinear Results... 170 59. 13.89 Hz Climatology Results.. 171 60. 13.89 Hz Transmission Loss.... 172 61. 111.1 Hz Bilinear Results ....173 62. 111.1 Hz Climatology Results.. 174 63. 111.1 Hz Transmission Loss 175 12 I- INTRODUCTION The growth of ocean acoustics as a theoretical and experimental science can be traced to a fundamental change in the conduct of warfare. During World War II, the belligerents employed scientific technology for the conduct of war to a degree theretofore unknown. In the United States scientific research was mobilized on a scale comparable to the mobilization of manpower and industry. Whole fields of scientific research were expanded to produce a technology for the conduct of the war. A portion of this research, the study of propagation of sound in the sea, was a response to a particularly effective (and technical) enemy weapon - the U-boat. The first wartime studies of underwater sound propagation hinted at sound transmission over long ranges within a natural oceanic waveguide known as the SOFAR channel. A series of experiments was conducted in which explosive signals were propagated in the SOFAfi channel over distances of up to 900 nautical miles (nm) [ 1 ]. In subseguent experiments explosive signals were propagated over paths of up to 2,300 nm [2]- Since World War II several uses for SOFAF. propagation have been devised, including the location of missile impacts within wide oceanic areas and the study of the low-f reguency ambient noise field. To predict and aid in understanding SOFAR propagation, a number of theoretical models for the transmission of sound in the ocean have been developed. The theoretical treatments have fallen generally into two broad catagories; ray-tracing and normal mode techniques. Ray-tracing techniques are based on the EiJconal equation, 13 2 2 VA • VA = N = [c /c (z) ] , (1.1) o where A is the ray function, N an index of refraction, c is o a reference sound speed, and c(z) sound speed as a function of depth. The Eikonal equation is restricted to high frequencies in which the sound speed fluctuation is small over one wavelength.. Normal mode techniques, more suited to low frequency sound, ace based upon the depth-dependent flelmholtz equation, 2 2 2 [V + u /c(z) ] * = 0, (1.2) where ¥ is the acoustic velocity potential and u is the angular frequency of the acoustic disturbance. A. DEVELOPMENT OF THE NORMAL MODE TECHNIQUE Since the work to be presented is based upon the normal mode technique, a brief overview of the development of underwater sound transmission models based upon normal mode theory is in order. 1 . Early Theory During World War II Pekeris developed two mathematical approaches which became the foundation for much of the later work. In 1946 Pekeris published a treatment of sound propagation in a layer characterized by a variety of sound speed gradients [3]. Although he was primarily concerned with sound penetration into the shadow zone caused by a negative sound speed gradient, he devleoped asymptotic solutions to the wave equation for a number of sound speed gradients. For the case of a linear sound speed gradient, c = c (1 ♦ az) , (1-3) o 14 he derived a solution which reduced to a Fourier- Bessel Series, which converged rapidly in the shadow zone and slowly in the insonified zone. In addition, he investigated a number of variable sound velocity profiles which included, 2 2 -1 c = c (1 - az) , (1.4) o for which the solution is a combination of Hankel functions of order one-thitd. In 1948 Pekeris published a method of determining the normal modes and resultant sound field in both two-and-three-layer isovelocity media [4]. To evaluate the sound field he applied a Hankel transform in the complex horizontal wave number plane to a vertical wave function U, which resulted in an expression for the acoustic velocity potential Y = fuj k dk . (1.5) (r o r r The wave function displayed a number of singularities and a branch cut. The solution consisted of a summation of normal modes, represented by a series of simple poles, and a ground wave in the lowest layer, represented by a branch cut (see Fig 1) . The branch cut extended from the branch point corresponding to c , the sound speed in the bottom, to i b This particular cut necessitated the evaluation of additional complex- valued poles between the branch cut and the positive imaginary axis (the so called "leaky modes") . The "leaky" modes were seen to be important only at short ranges. Pekeris evaluated the integral by the calculus of residues for the poles and the evaluation of a branch cut associated with the sound speed in the lowest layer. There were two significant modifications to the Pekeris technique which soon followed. Tolstoy reformulated the characteristic equation for the modal eigenvalues 15 Complex kr Plane - ^>*jWMvvo--VgMgfigA8fcAsW r • Real poles © Complex poles (Pe ken's) waxxo Branch cut (Pekeris) AVWo Branch cut (Ewing, et a I ) FIGURE 1 - Path of Bessel Integrals in Complex Plane 16 (location of the poles in the complex wave number plane) in terms of reflection coefficients for both upgoing and downgoing waves [5]. Tolstoy considered the effects of horizontally stratified layers above and below a depth of interest, z . A wave travelling up past z will be o o reflected with a reflection coefficient Rf. Similarly a plane wave travelling downward past z will have a o reflection coefficient S4- . Tolstoy showed that for an undamped mode representing the propagation of energy in the waveguide, the reflection coefficient above and below z must be related by Rf Rf = 1. (1.6) Rf and R^ are functions of the sound speed profile and a horizontal parameter (such as the horizontal wave number or the horizontal phase speed). Since Eg. (1.6) is satisfied only for certain discrete values of this parameter, it represents a characteristic equation for the modal eigenvalues. This general characteristic equation enabled normal mode analysis to be applied to a wide variety of environmental conditions. The similarity between the normal mode problem in acoustics and the SchrOdinger equation in quantum mechanics was recognized and provided a tool for further development. Bremmer published an analysis of the Wentzel-Kramer- Brillioun (WKB) method of quantum mechamics as an approximation to wave propagation in an inhomogeneous medium, opening the way for the use of the WKB solutions and characteristic equations as solutions to the vertical wave equations [ 6 ]• The second change to the Pekeris technique came when Ewing, Jardetzky and Press [7] suggested an alternate path for the branch cut (Fig 1). The Ewing, Jardetzky and Press 17 branch cut extended from the same branch point, to the origin, and then along the imaginary axis. This eliminated the need for evaluation of the complex-valued poles, as they were shown to be displaced to an adjacent Riemann sheet. The Pekeris model and its descendents were used in many theoretical applications, including shallow water propagation over fluid and elastic bottoms [8 - 15]; and was also used to study the effects of attenuating [16], sloping [17], and rough bottoms [18] on mode propagation. Up to this point a major limitation for normal mode and many ray-tracing techniques was the required assumption that sound speed be a function oniy_ of depth (to separate the range and depth-dependent equations) . In the ocean the sound speed profile varies gradually (and at times abruptly) with range. The applicability of the normal mode technique and, to a lesser extent, the ray-tracing technique over a horizontally varying medium was broadened by Weston in 1958 with the development of mode and ray invariants [19]. Weston observed that a mode may be assumed to keep its identity for slow horizontal changes in the sound speed profile. This is known as the adiabatic assumption (after its analogue in time-dependent perturbation theory of quantum mechanics) , implying that the energy propagated by a particular mode remains associated with that mode. A more recent paper by Milder expands upon the concepts of mode and ray invariants and derives a test for the applicability of the adiabatic assumption [20], The adiabatic assumption allows local horizontal stratification in an environment which changes slowly in the horizontal. Thus the normal mode method now became applicable to a wider variety of situations provided the horizontal oceanic variation was gradual over the long ranges. In addition, although not strictly required, the assumption of local horizontal stratification has been found useful in some ray-tracing work [21 ]. 18 2. Numerical Methods In the late 1950*s and early 1960fs there began to appear a number of numerical normal mode programs which took advantage of the digital computer. Many of the models were based at least in part on known solutions for certain classes of sound speed profiles. In 1959 Tolstoy and May developed a normal mode computer model based upon a series of horizontal layers in which the inverse of the square of the sound speed varied linearly with depth 2 2-1 c = c (1-az) . (1.7) o In such layers Tolstoy and May approximated the mode solutions (Bessel functions of order one-third) with twelfth-order polynomials [ 22 ]. For a given horizontal wave number they started with both a surface and bottom boundary condition, then iterated from the boundaries to some interior depth. The Bessel functions were matched between layers by invoking continuity of acoustic pressure and vertical particle velocity. At the interior depth the slopes of the solutions iterated from the surface and from the bottom were compared. Only certain values of the horizontal wave number, which defined the eigenvalues, gave first order continuity across the interior depth. The solutions thus formed defined the normal modes. Peter Hirsch studied axial sound focusing in the SOFAB channel using an analytic representation for the sound speed profile [23]. He assumed a sound speed profile -2 parabolic in c and symetric about the channel axis, 2 2 2 2-1 c = c (1-a z ) . (1.8) o For a pulse source this representation of the SOFAR channel gave successive focusing at about 6 nm intervals. Williams 19 and Home did another study of SOFAR channel axial focusing using a polynomial representation for the sound speed profile [24 ], 2 2 2 3 4-1 c = c (1-az-bz -cz -dz ) . (1-9) o This polynomial representation of the sound speed profile gave a better fit to the actual SOFAS channel, allowing representation of its asymmetrical character. Solving for the eigenvalues and normal modes, they described the sound field for a continuous wave (CW) source. Williams and Home found a much smaller degree of axial focusing than Hirsch; they also determined that focusing effects were sensitive to the detailed structure of the sound speed profile. Deavenport studied axial focusing using an Epstein profile, -2 2 c = A sech (z/h) + B tanh (z/h) + D, (1-10) with solutions which are hypergeometric functions [25]. He found a greater degree of focusing than Williams and Home, but less than from Hirsch* s symmetric profile. Bucker and Morris studied sound transmission in the surface duct, also using an Epstein profile for the sound speed [26]. The results compared favorably with transmission characteristics derived from experiment and ray-tracing models. In 1970 Williams published an outline of a normal mode finite difference technique [27] which started from either one or both boundaries. Using a finite difference equation based upon the vertical Helmholtz equation, the method iterated across the water column to either the other boundary or an internal depth. An iterative approach was used to find the mode eigenvalues by satisfying either the opposite boundary condition, or continuity of acoustic pressure and vertical velocity at the internal depth. A number of other models used this technique, - notably the 20 model of Newman and Ingenito [28], and the model by Kanabis [29]. The Newman and Ingenito model started with a fluid boundary condition at the bottom and used a finite difference technique to compute the eigenf unction across the water column to the surface. The eigenvalues were searched by finding those values of the horizontal wave number which most nearly produced an eigenf unction of zero at the surface. The Kanabis model was similar except that it used a higher order finite difference equation and employed the adiabatic assumption for horizontal variation. 3- Current Models In May 1973 the Acoustic Environmental Support Detatchment (AESD) of the Office of Naval Research held a workshop on non-ray-tracing acoustic propagation techniques [30]. The results from seventeen acoustic models were compared for calculations using a variety of sound speed profiles. A synopsis of the models and their results is given by Ref. [30]. A distinction was made between purely normal mode (in the sense of a complete set of orthonormal functions) and residue series models (where the modes are derived from the poles of the wave function) . At the AESD workshop a new model which represents a distinct departure from previous methods was introduced. The parabolic equation method, outlined by Fock for electro-magnetic propagation in the atmosphere [31], assumes the wave function to be a space dependent function, V, multiplied by a Hankel function. The assumed solution is applied to the Helmholtz equation, and terms in the second derivative of 7 with respect to range (r) are dropped (this assumes that the waves are predominantly radial - that is the corresponding rays are at relatively shallow angles) . The following parabolic differential equation results: 2 idV * J_d^v + JsTN (r,z)-1]V = 0, (1.11) oT 2ko*z* T 21 where k is an average wave number and N is an index of refraction. Since this is a parabolic differential equation, with a given starting solution, and surface and bottom boundary conditions, the solution can be solved in successive range steps. The method allows full range dependence for the sound speed profile. B. OPERATIONAL MODELS The U. S. Navy has a need to calculate low-frequency propagation loss over long ranges (in excess of 1000 km) . The Navy currently has variations of four standard operational transmission loss models and a number of experimental models. The operational models include the following: ■ FACT - The Fast Asymptotic Coherent Transmission Model uses modified ray theory and is primarily for passive sensor applications at low to intermediate frequencies (50 to 2000 Hz) . The FACT model computes propagation loss out to ranges representing several convergence zones [32]. FACT has been extended in an experimental version with multiple profiles and a bathymetry composed of linear segments along the line of bearing. The ray invariants of D. M. Milder [20] are used to connect rays between profiles. a NISSIM - The Navy Interim Standard Surface Ship Model is a ray theory model which is used primarily for short ranges and active sensors [33,34]. It has had some limited application to lower frequency passive sensors. ■ RP-70 - This model is a multiple profile ray-tracing model which uses linear sound speed gradients in both the vertical and horizontal directions. It is used for passive sensors and great ranges. ■ APE - The Acoustic Parabolic Equation model, based 22 on the parabolic aquation method (Eq. (1.10)), has recently become an operational model. It is primarily used for very low frequencies and highly range dependent environments. C. SCOPE OF THIS WORK 1 • Objectives The objective of this work is the creation of a computer model to estimate the transmission loss of low frequency CW sound over long ranges and through varying sound speed profiles. The model is intended to be relatively fast, sacrificing exactness (in the mathematical sense) for computational speed. The underlying theory is based heavily upon the WKB (Wentzel-Kramer-Brillioun) approximation with the eigenvalues, eigenf unctions, and normalization all based upon the WKB formulae. The model is intended for the following tasks: a Calculation of the detailed (fully coherent) transmission loss for ranges up to a few hundred nautical miles. In this case any sound speed variation with range must be slow enough to permit use of the adiabatic assumption. ■ Calculation of the incoherent or averaged transmission loss for long ranges (thousands of nautical miles) . 2. Departure from previous Works The following are a number of the major differences between this work and previous normal mode models: ■ The sound speed profile is assumed to consist of layers (up to 50) in which the sound speed varies linearly 23 2 in c (Eq- (1.6)). Previous models have used the analytic solution for such layers to form the eigenf unctions (modes) and find the eigenvalues (Tolstoy and May [22], Gordon and Petersen [30], and Bartberger and Ackler [35]). 2 In this model the c linear profile is primarily used for the ease with which the vertical wave number may be integrated - an operation required for many of the WKB calculations. The analytic solution (the Airy functions) is used as a mode solution only where the WKB approximation breaks down. ■ The sound speed profile is divided into separate ducts corresponding to different families of modes. The calculation of eigenvalues and modes is carried out family by family. ■ Within each family the behavior of the eigenvalue as a function of mode number is approximated by a fifth-order polynomial. This polynomial provides the first estimate of the eigenvalues for each mode. For calculation of the averaged or incoherent transmission loss at great range, this estimate is sufficiently accurate for computation of the eigenf unctions. For detailed transmission loss at closer range, the eigenvalue is refined by a maximum of five iterations. ■ The WKB approximation has been extended to incorporate the effects of a . barrier between multiple sound channels, and the effects of the surface and bottom boundaries. The eigenf unction is found by a combination of WKB, Airy and extrapolated functions. ■ A range dependent environment has been included by the use of the adiabatic assumption. A previous difficulty with the adiabatic assumption for normal mode use has been the ambiguities for modes in different sound channels of multiple channel profiles. This difficulty has been overcome with the use of a formalism which 24 matches modes within equivalent channels and approximates the transition of modes between channels. D. REFERENCES 1 Ewing, Maurice and J. Lamar Worzel, "Long Range Sound Transmission", in Propagation of Sound in the Ocean, Geol. Soc. Am. Memoir~277~l5~0c£ T^UB. " ~ 2 Herring, C, "Transmission of Explosive Sound in the Sea", in Physics of Sound in the Sea, National Defense Research Council, Summary TecEnical Report of Division 6, Vol. 8, 1946. Reprinted by Headquarters Naval Material command, Washington, D.C., 1969. 3 Pekeris, C. L., "Theory of Propagation of Sound in a Half-Space of Variable Velocity under Conditions of Formation of a Shadow Zone", J. Acoust. Soc. Am., Vol. 18, No. 2, Oct 1946, pp. 295-315. 4 Pekeris, C. L., "Theory of Propagation of Explosive Sound in Shallow Mater", in Propagation of Sound in the Ocean, Geol. Soc. Am. Memoir~27, T5~0"cT T9"48\ " 5 Tolstoy, Ivan, "Note on the Propagation of Normal Modes in Inhomogeneous Media", J. Acoust. Soc. Am., Vol. 27, p. 274, Mar 1955. 6 Bremmer, H. , "The W.K.B. Approximation as the First Term of a Geometric-Optical Series", Commun. Pure Appl. Math., Vol. 4, 1951, pp. 105-115. 8 Tolstoy, I. , "Dispersion and Simple Harmonic Point Sources in Wave Ducts", J. Acoust. Soc. Am., Vol. 27, p. 8 97, Sep 19 55. 7 Ewing, W. Maurice, Wenceslas S. Jardetzky and Frank Press, Elastic Way.§§ in Layered Media, McGraw-Hill, 9 Tolstoy, I. , "Shallow Water Propagation Tables, Three Layer Waveguides", Hudson Labs, Columbia Univ. , Tech Rept 75, Dec 1960, AD-256-699. 10 Clay, C. S. , "Propagation of Band Limited Noise in a Layered Waveguide", J. Acoust. Soc. Am., Vol. 31, p. 1473, Nov. 1959. 11 Tolstoy, I.. "Guided Waves in a Fluid with a Continuously Variable Velocity Overlying an Elastic Solid: Theory and Experiment", J. Acoust. Soc. Am. , Vol. 32, p. 81, Jan 1960. 12 Williams, A. 0., Jr., "Some Effects of Velocity Structure on Low-Frequency Propagation in Shallow Water", J. Acoust. Soc. Am., Vol. 32, p. 363, Mar 1960. 13 Williams^ A- 0., Jr., and R. K. Eby, "Acoustic Attenuation in a Liquid Layer over a *Slow' Viscoelastic Solid", J. Acoust. Soc. Am., Vol. 34, p. 836, June 1962. 25 14 Backer/ H. P. , "Normal Mode Sound Propagation in Shallow Water". J. Acoust. Soc. Am. , Vol. 36, p. 251, Feb 1964. 15 Bucker, H. P., and H. E. Morris, "Normal-Mode Intensity Calculations for a Constant-Depth Shallow Water Channel". J. Acoust. Soc. Am., Vol. 38, p. 1010, Dec 1965. 16 Kornhauser, E. T., and W. P. Raney, "Attenuation in Shallow Water Propagation Due to an Absorbing Bottom", J. Acoust. Soc. Am., Vol. 27, p. 689, July 1955. 17 Eby, R. K., A. 0. Williams, Jr., R. P. Ryan, and Aul Tamarkin, "Study of Acoustic Propagation in a Two-Layered Model", J. Acoust. Soc. Am., vol. 32, p. 88, Jan 1960. 18 Clay- C. S., "Effect of a Slightly Irregular Boundary on the Coherence of Waveguide Propagation", J. Acoust. Soc. Am., Vol. 36, p. 833, May 1964. 19 Weston, P. E., "Guided Propagation in a Slowly Varying Medium", Proc. Phys. Soc. London, Vol 73, p. 365, 1958. 20 Milder. D. Michael, "Ray and Wave Invariants for SOFAR Channel Propagation", J. Acoust. Soc. Am., Vol. 46, p. 1259, 1969. 21 Smith, P. W.- "Averaged Sound Transmission in Range-Dependent Channels", J. Acoust. Soc. Am., Vol. 55, p. 1197, June 1974. 22 Tolstoy, I. and J. May, "A Numerical Solution for the Problem of Long Range Sound Propagation in Continuously Stratified Media, with Applications to the Deep Ocean", J. Acoust. Soc. Am., Vol. 32, p. 655, June 1960. 23 Hirsch, Peter, "Acoustic Field of a Pulsed Source in the Underwater Sound Channel", J. Acoust. Soc. Am. , Vol. 38, p. 1018, Dec 1965. 24 Williams, A. 0., Jr., and William Home, "Axial Focusing of Sound in the SOFAR Channel", J. Acoust. Soc. Am., Vol. 41, p. 189, June 1967. 25 Deavenport, R. L. , "A Normal Mode Theory of an Underwater Acoustic Duct oy Means of Green's Function", Radio Science, Vol. 1, p. 709, 1966. 26 Bucker, H. P., and H. E. Morris, "Epstein Normal-Mode Model of a Surface Duct", J. Acoust. Soc. Am., Vol. 41, p. 1475, June 1967. 27 Williams, A. 0., Jr., "Normal Mode Methods in the Propagation of underwater Sound", in Underwater Acoustics, ed. by R. W. B. STepHens, T7iley-TnTer science, 1970. 28 Newman, A. V.. and F. Ingenito, A Normal Mode Computer Program for Calculating Sound Propagation in Shallow ¥aTer wiTTS an Ir]5trary Velocity Profile, ^TTavaT ResearcH Xa5 HemorancTum^BeporT 23"HT,"~Jan T97Z. 29 Kanabis, W. G. , Computer Programs to Calculate Normal Mode Propagation and~Appricarions~€o tEe Xnalysis of Explosive souncT in "the BTTT Ianae,"~"NavaI Underwater Sys'Eems^en'Eer TecEnicaI~Repor£ 53T97 6 Nov 1972. 26 30 Spofford, C. H.r A Synopsis of the A ESP Workshop on Acoustic Propagation Mojjglllng. Using ""Nqn'-gay-Tracing Tecanigues, "Scoustic E"nvironmentaT'~5"uppor"E~"Detatcniiient ISESUr^Technical Note TN-73-05, Nov 1973. 31 Focic, V. A., EI§ctromagnetic Diffraction and propagation Problems, Pergamon Press, 1965, pp. 2T 3- 235". 32 Spofford, C. «., The FACT Mo.de.1. Vol. I, AESD, Maury Center for Ocean Science, ffC Heport 109, Nov 1974. 33 DELETED. 34 Watson, W. H., and R. W. McGirr, An Active Sonar ps£JL2r5L§Hce Prediction Model, Naval Undersea "Researcn" and UeveTopmenF" CenEer ~7flUC) Technical Publication NUC-TT-286, April 1972. 35 Eartberger, C. and L. Ackler, Normal Mode Solutions and Computer Programs for Underwater SounH Pf opaqar4on, Part T, Uaval" Air"" Development Center Report no. NADC-72001-AE, 4 April 1973. 27 II. TH^OaYi HOMOGENEOUS WAVE EOJJAT.JON In what follows it is assumed that the ocean is a horizontally stratified medium of constant density, with a perfectly flat and pressure release surface, and a flat bottom. The input data for most ocean acoustic models represent a linearly-segmented sound speed profile. The linear segments do not represent the detailed structure of the ocean sound speed profile; they are a convenience for interpolating between the data. In ray trace models this linear gradient takes on added significance because of the simplified ray path thus formed. Because of the analytic solutions available, when using normal mode techniques it is convenient to approximate the sound speed structure between -2 data points as linear in c . The profile will be assumed to consist of layers within which the sound speed varies -1/2 c(z) = c [1 - a (z-z ) , z < z < z , (2.1) i i i i+1 where c and z are the sound speed and depth at the top of i i th the i layer (Fig 2) . If u is the angular frequency of the acoustic disturbance, then the sound speed profile can be represented in terms of the wave number, 3c. For the profile 2 of Eg. (2.1) k varies linearly with depth, and the gradient 2 of k , y , is defined by 2 2 2 2 jc = a) /c = k -y (z-z ) , z < z < z . (2.2) i i i i+1 28 SOUND SPEED FIGUHE 2 - Sound Speed Profile 29 The horizontal wave number k and the vertical wave number r k are related by z 2 2 2 k = k - k . (2.3) z r The vertical wave number is real for k > k , and is r imaginary for for k < k ; r k =-i8, z where 2 2 2 6 = k -k (2.4) r for k < k . r The theoretical development which follows can be divided into five steps outlined below: ■ The homogeneous wave eguation is established, and from it is separated a depth-dependent equation. The depth-dependent eguation must be solved in terms of eigenvalues and corresponding eigenf unctions (the normal modes) . ■ Two approximate solutions to the homogeneous depth-dependent equation, the HKB solutions, are derived. One oscillatory solution corresponds to the trapped field in the waveguide, the other, exponentially decaying, represents the diffraction of sound at the edge of the waveguide. The HKB solutions fail where they approach each other at the edge of the channel (turning points) . From one of the 3KB solutions a characteristic equation is developed, whose solutions are the modal eigenvalues. ■ A third approximate solution to the homogeneous equation is derived. This solution is used to connect the 30 two WKB solutions. As a consequence, for each eigenvalue an unnormalized eigenf unction can be determined as a function of depth. * The treatment of the homogeneous wave equation is completed with consideration of the effects of the surface and bottom boundaries. ■ The final step is the introduction of an inhomogeneous term (representing a point source) into the wave equation. The inhomogeneous equation is then solved by the residue series of a Green's function which is formed from the previously developed homogeneous solutions. Solution of the inhomogeneous equation allows determination of the mode amplitude to match a point source. A. HOMOGENEOUS WAVE EQUATION The homogeneous lossless wave equation for velocity potential ¥ is V2Y - 1 d£ I = 0, (2.5) c* eft* The acoustic pressure, P, is related to the velocity potential by P = - p d¥ . For a monofrequency acoustical disturbance with time dependence exp (-i cot), P = i copy. The boundary condition at the surface requires that Y (z=0) = 0, (2.6) where z is the depth below the surface. The boundary 31 condition at the bottom z can be written in the form b 1 AX = Jii {2.7) y ciz where y will be determined by the particular bottom model used. In addition, certain restrictions must be placed on 1. ¥ must represent an outgoing wave in the horizontal direction and an exponentially decaying wave at great depth, so that |^| -* 0 as r + °° or z -* °° . 2. There must be no discontinuity in pressure and particle velocity in the water column so that p * and the gradient of * must be continuous. 1 • The Helmho It z Eg.ua t ion If the time dependence of y is explicitly stated, V = $ exp (-iujt) , (2.8) where $ is a function only of space, then substitution of Eg. (2.8) into Eg. (2.5) yields V2 $ ♦ w* = 0. (2.9) Expanded in cylindrical coordinates with the origin at the sea surface and z positive downward, Eg. (2.9) becomes 1 d_(rd£) + d£$ + us =0, (2.10) r clr clr 3z? c? for cylindrically symmetric motion. 2. The Hankel Transform One way to separate the depth dependent terms in Eg. (2.10) is to apply a Hankel transform. The Hankel transform variable, k , represents a horizontal wave number, r 32 and transforms the solution from a depth and range space (z,r) to a solution in depth and horizontal wave number space (z,k ) . The transform pair is defined r U(z,k ) = U (z,r) J (k r) r dr, (2.11) r qJ or and (z,r) = / C(z,k ) J (k r) k dk . (2.12) Qj r o r r r Application of this transform to Eq. (2.10) yields 2 d£0 + Im* " k ] U(z,k ) = 0. (2.13) az* c* r r This is a homogeneous differential equation with the following conditions imposed as a consequence of the conditions on y : a) U(0,k ) = 0, the free surface (2.14) r boundary condition. b) 1 dU = ]i, at the bottom, where z = z . U 3z b c) As z->-°° , 0 must represent a decaying wave. d) p U is continuous. e) dU is continuous. 3z 3. Nature of the Solution The function U(z,k ) defines one normal mode, which r represents the depth component of an acoustic wave propagating outward with cylindrical divergence. The horizontal wave number defines the phase speed, c = oj / k . (2. 15) P r The phase speed represents the speed with which surfaces of constant phase are propagated in range. 33 Inspection of Eq. (2.13) reveals some qualitative properties of the solution. For k > k the vertical wave r number, k is real and the solution is expected to be z oscillatory in nature. For k < k the vertical wave number r is imaginary, and the solution is expected to be exponential in nature. Regions for which k is real (cc ) form an z p acoustic channel or waveguide. Within such a channel the acoustic waves are refracted between depths at which k =0, z or reflected from boundaries. The ocean normally forms an acoustic waveguide for phase speeds less than the speed of sound in the bottom, c . b At these phase speeds the modes correspond to waves which are either completely refracted or strike the bottom at a grazing angle less than the critical angle. For the corresponding values of k (k >w/c =k ), Eg. (2.13) with the r r b b conditions of (2. 14) forms a Sturm-Liouville boundary value problem over the interval z=0toz=z. As a consequence b there is a finite number of eigenvalues k for the r,n parameter k and corresponding eigenf unctions, U , which r n solve the Sturm-Liouville problem [ 1 ]. Thus the range of k r > k is termed the discrete spectrum, b At phase speeds greater than the sound speed in the bottom, the differential equation and boundary conditions are met for all phase speeds or k . For these high phase r speeds there is a continuous and infinite set of solutions 34 for all k < k . r b The reader familiar with the SchrBedinger equation of quantum mechanics will recognize that the sound velocity profile is an analogue to a potential well. The sound speed profile itself corresponds to the potential function, the phase speeds correspond to the energy levels, the surface boundary to an infinitely high "wall", and the bottom boundary to a "wall1* with height that corresponds to the sound speed in the bottom [2] [3]. B. THE WKB METHOD In this section the WKB method, an approximate solution originally developed in quantum mechanics, is applied to solution of the wave equation for normal modes in the ocean. 1. The WKB Solution To derive the WKB formulae a method used by Schiff is adapted to oceanic parameters [4]. By defining a scale parameter, h = 1/w , and an index of refraction, 2 2 2 N = [1/c - 1/c ], P the vertical wave number can be re-written, 2 2 2 k = N / h . (2.16) z Assume that the depth function has the form U = A exp[iS(z}_], (2.17) 35 where S is some (Unknown) function of depth. Substitution of Eqs. (2.16) and (2.17) into Eg. (2.13), gives 2 2 ihS» ■ - S» + N = 0. (2. 18) Expand S in powers of h, S = S + hS ♦ ... ; (2.19) 0 1 and equate terms in the first two powers of h, resulting in two equations 2 2 -5« + N = 0, (2.20) and iSfl - 2S»S» = 0. 0 0 1 Integration of the first equation yields f z J N dz'r S (z) = ± / N dz'r (2.21) Z0 and substitution into the second gives S (z) = i log N. (2.22) 1 1 e Substitution into Eg. (2.17) gives the two WKB solutions. For k > k , r -1/2 a a 0(z) = k {a cos( k dz) + b sin ( k dz) } . (2.23) z o 7J z o 7 J z z0 Z0 This solution will be referred to as the HKB oscillatory solution and can also be written -1/2 rz 0 (z) = k cos( k dz-0 ), (2.24) I z J z o o where 36 6 = arctan (b /a ) . (2.25) o o o For k < k the HKB exponential solution is obtained, r z z -V2 r r 0 = {c exp( B dz) + d exp(- 6 dz) } . (2.26) II o J o J zo 2. ffKB Validity Requirements In order to expand S in powers of h (Sq. (2.19)) successfully and to separate the first two terms, the second term, hS , must be small in comparison to the first, S . 1 0 This requirement is met if [ 5 ] -2 IhSJ/SM = Ilk d_ k I « 1. (2.27) 10 Z z 3z z A WKB validity factor, F , will be defined as equal to the w two-thirds power of the expression in Eq. (2.27). (The two-thirds power is used to facilitate comparison with the Airy validity factor derived in a later section.) From Eq. (2.27) with the help of Eq. (2.2) 2/3 -1 -1/3 F = 1 |z-z | v << 1, (2.28) w IT t where z is the turning point depth. Thus |F | must be much t w less than one for the WKB solution to be a good representation. In terms of the horizontal wave number the expression for F becomes, w 2/3 2 2-1 2/3 F = 1 |k -k j v • (2.29) w tf r 37 3- Properties of the WKB Solutions The solution to Eg. (2.13) is characterized by an oscillatory function (Eg. (2.24)) at those depths for which k is real. This function corresponds to the summation of z an upgoing and a downgoing vertical wave. At those depths for which k is imaginary, the solution represents the z lateral wave in a shadow zone. For the reference depth z o in Eg. (2.24) and (2.26), the turning points or a reflective boundary will be used. The Q of Eg. (2.24) represents a phase change as o the vertical wave is refracted at the turning point or reflected at a boundary. Considering the oscillatory solution as the combination of an upgoing and a downgoing wave, 9 represents the phase change for only one of the o waves. Thus the total phase change experienced by a wave as it travels through a turning point or is reflected from a boundary is twice 9 . o At the turning point F becomes singular; thus there w is some region bounding the turning points for which the WKB solutions are not valid. 4. WKB Characteristic Equation The homogeneous differential eguation (Eg. (2.13)) can be solved and both boundary conditions satisfied only for certain eigenvalues of k . After developing the WKB r approximations, the next step is to use the WKB oscillatory solution to determine a characteristic eguation for the mode eigenvalues, k r#n 38 Consider some arbitrary depth, z , between the two 1 turning points for a specific k (Fig 3) . Assume r,n U(z,k ) is a solution to Eg. (2.13) which meets all the r,n conditions of (2.14). From Eg. (2.24), at depths above z the solution must be -1/2 r1 0 (z -) = k cos ( k dz-e ) , 11 z J z u ztu where z is the upper turning point depth and 6 is the WKB tu u phase angle at the upper turning point. Conversely, if the lower turning point z is used as the reference level, the tl solution below z must be 1 Ztl -1/2 U (z ♦) = k cos ( 2 1 z zi For the two expressions to be equivalent, they must be k dz-9 ) . J z 1 continuous through the first derivative at z . This 1 requirement is met only if / zti k dz-Q -9 = n tt , n=0,1,2, (2.30) z u 1 ztu This is the characteristic equation for the WKB approximation. If d and 9 equal tt /4 this characteristic u 1 equation becomes the 3ohr-Soramerf eld quantization rule for a potential well [6]. Consider a wave which propagates upward from z to the upper turning point z , is refracted downward to z , tu tl and then is refracted upward back to depth z . The 1 39 SURFACE SOUND SPEED WKB oscillatory WKB exponential upper turning point lower turning point WKB exponential FIGURE 3 - The WKB Solutions 40 characteristic equation implies this wave has undergone a phase change which is an exact multiple of 2 C. REFERENCES 1 Churchill, Ruel V. , Fourier Series and Boundary Value Problems, (2nd. edition)*, HcSraw-Hill, 19^3, pp. 67=757"" 2 Schiff, L. I.; Quantum Mechanics, McGraw-Hill, 1955; pp. 27-34. 3 Williams, A. 0., Jr.; Normal Mode Methods in Propagation of Underwater Sound; Underwater Acoustics, ed. by R.W.B. Stephens, Wiley Interscience, pp.337 4 Schiff, L. I.; ibid., pp. 184-186. 5 Schiff, L. I.; ibid., pp. 186-187. 6 Tolstoy, I. and C. S. Clay, Ocean Acoustics, McGraw-Hill, 1966, pp.51. 41 III. THEORY: AIRY SOLUTION AND CONNECTION FOBMULAE After deriving tha WKB solutions, which are good approximations away from the turning points, the next step is to derive a solution valid at the turning points. This solution, called the Airy solution, will be used to match coefficients between the two sets of WKB eigenfunctions determined from Eg. (2.23) and (2.26). The derivation will follow the method used by Lauvstad [1]. A. THE AIBY SOLUTION 2 For some region about the turning point, k varies r linearly with depth and Eg. (2.13) may be expressed d£U - y (z-z )0 = 0, (3.1) dzz t where y is defined in Eg. (2.2) and z is the turning points With the introduction of a new variable, 1/3 x = Y (z-z ), (3.2) the second derivative with respect to depth becomes 2/3 d£ = d (d_dx) = y <*£ » dz? cfz Uxctz 3x7 resulting in the Airy eguation, d2U - xU = 0. (3.3) Hx* H2 Solutions to this differential equation are Bessel functions of order one-third. The solution to Eg. (3.3) may be expressed in terras of a particular combination of these Bessel functions, the Airy functions Ai and Bi [ 2 ] U(x) = A Ai(x) + B Bi(x). (3.4) o o With approximate expressions for the Airy functions for large and small argument, the two WKB solutions may be matched to the Airy functions and, in the process, matched across the turning point. 1 • Airy Solution for Small Argument For small x the expressions for Ai and Bi can be represented as ascending series [3]- Taking the first two terms of those series, we have for |x| << 1, U (x) = C (A */lB ) ♦ C (/JB -A ) x, (3.5) III 1 o o 2 o o where c = 0.35502805, c = 0.25881940 . 1 2 This expression, the Airy solution, is a good approximation at the turning point, precisely where the WKB solutions break down. 2- Airy Validity Factor The WKB solutions are good approximations at some depth beyond the turning points. In addition, there is some interval about the turning point z within which the Airy solution is a good approximation. The requirement that |x| << 1 in Eg. (3.5) may ba imposed through the use of an Airy validity factor, F : a 43 1/3 F = |x| = J y (z-z ) | « 1, (3.6) a t or -2/3 2 2 F=|xj=Y Ik - k |. (3.7) a r The validity factors F and F are plotted as functions of a a w V3 scaled depth, Y (z-z ) (Fig 4). If limits of 0.3 and 0.448 (0.3 to the two-thirds power) are used for F and F a w respectively, the Airy solution is a good approximation for x < 0.3 and the WKB solution is a good approximation for |x| > 0.87. There is a "gap" for 0.3 < |x| < 0.87 within which neither the Airy nor WKB solutions is a good approximation. ' 3. Airy Solution for Large Argument To determine the correspondence between the Airy solutions and a particular WKB solution, the coefficients of the Airy functions for large argument and the WKB solutions are matched. Introduce a variable (solely for notational convenience) , 3/2 ? = 2 x ; 3 for Jx| >> 1, the Airy solutions may be written in terms of the asymptotic expansions for Ai (x) and Bi (x) [4]. For the shadow zone region ( k > k) , the Airy solution, for r asymptotically large argument, is -1/4 -S C 0(x) = 1x ll A e + B e ]. (3.8) Vtt 2 44 °J PUD Md a* Q QJ H U to O' G o •H +J • u a 3 IV, w id •H U i •P •H •* H > W crj C!) H 45 For the oscillatory region ( k < k) , the Airy solution for r large argument is -1/4 U (-x) = 1x [B cos(£+it) + A sin (£ -ht ) ]. (3.9) JP o % o H B. CONNECTION FORMOLAE The next step is to equate both WKB solutions, U and U , to the Airy solution, U , then match the WKB II III solutions across the turning point. In addition, the case of a barrier (sound speed relative maximum between two channels) is considered, and the WKB solutions for the channel on one side of the barrier are matched with the solutions on the other side of the barrier. The WKB solutions are a function of the integral of the 2 vertical wave number. For a single linear segment (in k ) this integral is expressed z : |k |dz = 2 Y jz-z | = C . (3.10) 1/2 3/2 Ikjdz = 2 For notational convenience a constant, g, is defined by 1/6 -1/2 g = y ( * )• 1 . WKB Oscillatory and Airy Connection By matching coefficients, the oscillatory WKB solution (U ) and the Airy solution (U ) are related by I III 46 a = 1g (A +B ), b = 1g (A -B ), (3.11) o yz o o o /7 o o A = _1 (a +b ) , B = 1 (a -b ) . o /Zq 00 o 77g o o 2. WKB Exponential and A^ry, Connection The exponential HKB solution (D ) and the Airy solution (U ) are related by III d = q A , c = g B , (3.12) o 2 o o o A = 2 d , B=1c. o g o o g o 3. WKB Connection Formulae Finally, by combining Eg. (3.11) and (3.12), the two WKB solutions are matched by a = 1 (2d +c ) , b = 1 (2d -c ) , (3.13) o 77 00 o yT o o c = 1 (a -b ) , d = 1 (a +b ). o /7 00 o Z/T o o The % of Eg. (2.20) can now be re-defined as o 3 = arctan (b /a ) = arctan [ (2d -c )/(2d +c ) ]. (3.14) o 00 0000 C. BARRIER CONNECTION FORMULAE When a sound speed profile has two relative minima (two channels) , it is necessary to describe the lower turning point phase angle of the upper channel, 8 , in terms of the 47 upper turning point phase angle of the lower channel, 0 , and k (Fig 5). Consider a wave which approaches the lower r turning point of the upper channel, z . The wave is 1 refracted and vertexes in the vicinity of the turning point, z . Some of the wave energy "leaks" through the barrier and 1 enters through the upper turning point, z , into the lower channel. If the barrier is thick enough for the wave to be considered as progressing from an Airy representation to a WKB representation, and then back to an Airy representation, the connection formulae (Eg, (3.13)) may be invoked twice. Thus the WKB coefficients in the upper ensonified channel (a ,b ) in terms of the coefficients for the lower channel o o (a ,b ) can be expressed as 2 2 L -L a = (a -b )e + 1(a +b ) e , o 2 2 5 2 2 b = (a-b)e - 1(a +b )e , (3.15) o 2 2 5 2 2 where 7 2 2 ■/ Zl L = / 8 dz; [ 6 is the imaginary vertical wave number in the barrier, see Eg. (2.4) J. The MKB phase angle at the lower turning point of the upper ensonified channel, 9 (z ) , can be 1 1 expressed in terms of the phase angle at the upper turning point in the lower channel, Q (z ) : 2 2 43 SOUND SPEED \ "leaked" i Wave lower transmitted wave FIGORE 5 - Double Sound Channels and Barrier 49 d = arctan (b /a ) 1 o o L -L = arctan {[ (1-tanQ )e - t(l + tan-0 )e ] / L -L [ (1-tanS )e + 1 (1+tane )e ]} . (3.16) If the upper channel is considered the ensonified channel and the lower channel that into which energy is "leaked" , the relative amplitude of the lower solution represents the transmission coefficient for the barrier. This requires that the acoustic amplitude at the top of the barrier exceed the amplitude at the bottom of the barrier. Thus, the amplitude of the WKB exponential solution at the top of the barrier should be larger in magnitude than the WKB exponential solution at the bottom of the barrier. In terms of the coefficients of Eg. (2.26) this may expressed 2 L . -L 2 [c ♦ d ] >[ce +de } , o o o o or -L jc | < |d | e . (3.17) o o When this is substituted into the connection formulae (Eg. (3. 13) ) , a restriction is placed on the lower phase angle of the upper channel ■n - arctanjc /2d | < 0 < ir + arctan|c /2d |. (3.18) if o o 1 5 o o Thus for a wide barrier, the WKB phase angle on the ensonified side of the barrier approaches "f/H. As the barrier narrows, the phase angle may be perturbed an increasing amount on either side of the t/4 value. As the barrier vanishes, L in Eg. (3.17) approaches zero and the phase angle may vary from between 18.66° to 71.34°. 50 D. REFERENCES 3 a Lauvstad, V.R., An Expansion Technique for Rendering the HKBJ Method" CTnifofmlx 711137 ""CTATQ SACTYHTOSw BesearcI^CenTre- TecEnlcal ETeporfTJo. 203, 15 Dec 1971. Abramowitz, M. , and I- A. Stegun, Handbook of Mathematical Functions, National Bureau of STancIara's IppTied HatEematics" Series No. 55, U. S. Govt. Printing Office, 1964, p. 446. Abramowitz, M. and I. A. Stegun, ibid, p. 446. Abramowitz, M. and I. A. Stegun, ibid, p. 448. 51 IV. BOUNDARY CONDITIONS In this chapter expressions for the surface and bottom boundary conditions are developed in terms of the Q of o Eg. (2-24). Also, in order to evaluate the eigenf unction normalization factor (derived later) , expressions for the derivative of 0 with respect to k are obtained. The o r derivations will be divided into three cases, based upon whether the sound field at the boundary is best described by the WKB exponential solution (Eq. (2.26)), the Airy Solution (Eq. (3.5)), or the »KB oscillatory solution (Eq. (2.24)). A. THE SURFACE BOUNDARY The pressure release boundary condition at the surface as defined by Eq. (2.14a) is applied to one of the three solution cases. The coefficients thus derived are matched to produce a value for 6 in Eq. (2.30). u 1. Case I Consider a surface boundary which is in the shadow zone region. If Eq. (2.27) is satisfied at the surface, the WKB shadow zone solution (Eq. (2.26)) is a good representation at the surface (Fig 6) . Apply the surface boundary condition to obtain -V2 rtu rtu U (0) = 6 {c exp( gdz) + d exp(- Bdz) } = 0, (4.1) II o o J o J cr o 52 SOUND SPEED FIGURE 6 - Surface Boundary Condition - Case I 53 where g is the value of g at the surface. Using the 0 notation fztu .u = J 6 dz, (4.2) 0 and with the connection formulae (3.13} and (3.14), the surface boundary condition may be written, -2Lu -2Lu 0 = arctan (b /a ) = arctan [ (2+e )/(2-e ) ]. (4.3) u o o As the depth of the turning point increases, Lu increases in magnitude, and 9 approaches tt/4. This result corresponds u to the tt/4 phase shift for a refracted wave in an unbounded medium, as derived by many authors including Brekhovskikh £1], Schiff [2], and Tolstoy and Clay [3]. Differentiating with respect to k obtain r -2Lu -4Lu -1 d (3 ) = -8e [4+e ] d (Lu) . (4.4) die.- u ak„ 2. Case II If |x| << 1 in Eg. (3.5), and the surface and 2 turning point occur in the same linear segment of k , then the Airy solution (Eg. (3.5)) can be matched to the surface boundary condition (Fig 7) c (A +,/3B ) + c (/3B -A )x = 0. (4.5) 1 o o 2 o o Using the turning point connection formulae, the upper turning point phase shift may be expressed, 0 = arctan (b /a ) u o o 54 SURFACE SOUND SPEED FIGOEE 7 - Surface Boundary Condition - Case II 55 where b /a = [c (/3-1)x ♦ c (/T+1) ] / o o 2 1 [c (/3+1)x ♦ c (/3-1) ]. (4.6) For values of |x| < 0.3 the upper turning point phase shift increases from 65.31° to 87.00° as the turning point approaches the surface (and the mode becomes surface reflected) . This differs from the WKB exponential upper turning point phase shift, which varies from 45° to 72.57° as Lu decreases from a large number to zero. The derivative of 0 with respect to k is u r 2 2-1 d (G ) = (a +b ) 2k c [a (/3-1)-b (/S+1) ]- (4.7) cTTcu oo r2o o When x equals zero, which corresponds to a wave which is refracted with grazing incidence at the surface, 6 = arctan [ (/T + 1 ) / (/3- 1 ) ], u = 5tt/12 radians or 75°. This condition is analogous to the ray that grazes the surface at zero degrees, or is just tangent to the surface. This result can be duplicated by starting with the connection formulae given by Schiff [4], rather than with the connection formulae (3.13). Since Eg,. (4.5) is an approximation for U on both III sides of the turning point (for small values of x) , Eq. (4.6) should provide a good representation for small negative as well as small positive x. This corresponds to waves reflected from the surface at extremely small grazing 56 angles. 9 as a function of x increases as the denominator u of Eg. (4-6) approaches zero. This occurs for a value of x = [c (1-/31) ] / [c (1+/T) ] = -0.3675516 . 1 2 At this value, G = it/2 radians = 90°, u and we have the most common form of the pressure release boundary condition - a phase shift of n radians or 180 degrees (here split equally between the upgoing and downgoing waves) . 3. Case III In this case consider a wave that is reflected from the surface at an angle large enough so that the WKB representation is valid. Thus the boundary condition becomes (Fig 8) z -1/2 0 = k cos( / k dz-9 ) = 0. (4.8) / I z QJ z u At the surface, z = z and o 3 = tt/2 or 90°. (4.9) u This is the expected result for a pressure release surface. *• General Cases The results for the three cases described above are 2 shown in Fig 9 as a function of x, for k a linear function of depth. The HKB validity criterion F is also displayed w (the validity criterion for the Airy solution is |x|). The 57 I Q_ LJ O surface OUND SPEED PROFILE SOUND SPEED U, =0 -O s* WKB 1 Oscillatory hot ut ion I PHASE SPEED FIGURE 8 - Surface Boundary Condition - Case III 58 o CO O CO LU LU (T g * -i- f — ■♦ — i 6 o AIRY Solution WKB Solution t — i — | — i — i — i — i — | — i — i — r -0.5 0.0 X — J — I — I — » — I— 0.5 I.C -1,0 FIGURE 9 - Surface Boundary Condition - 9 vs. x u 59 Airy analysis shows that there are three distinct cases for the effect of the surface pressure release boundary upon the upper turning point phase shift. The extent of each case depends upon the gradient of k at the surface. Note that r there is a "gap" between those values of jxj which are sufficiently small for |F | << 1 , and those values which are a sufficiently large for |F | << 1. This gap could be w shortened by the inclusion of more terms in the ascending series for the Airy functions. In addition, the sound speed profiles normally do not allow the representation of the sound speed curve between the surface and the upper turning point as a single 2 linear gradient in k . Usually there is a combination of such linear segments, which may affect the results. Consider a profile near the surface as a combination of linear segments as in Fig 10. Since the boundary condition at the surface must be satisfied, it must be determined which, if any, form of the solution may apply at the surface. In most cases there should be a range of relatively low phase speeds for which the WKB shadow zone solution is a good approximation at the surface. This range corresponds to the deep channel refracted waves (Fig 10) . At the other end of the phase speed spectrum, there will be some range for which the WKB oscillatory solution (with a phase shift of u /2) is a good approximation. In addition there will be some narrow band of phase speed, centered about the surface sound speed, for which the Airy solution is a good approximation. The gaps between these ranges of good approximation could be bridged by using the full Airy function 60 surface x H Q_ a SOUND SPEED £> Uj CO o 1 to <© a o •H ■P •H (3 O U >i u «s £3 O CO QJ U 4-1 U CO H S33U93Q U! "8 63 As with the surface boundary condition the lower turning point phase change, 6 , will be investigated for the three forms of the solution. 1 . Case I For modes refracted at depths well above the bottom, the WKB exponential solution is a good approximation (Fig 12) . The derivative (with respect to depth) of the WKB exponential solution (Eg- (2.26)) is -1/2 LI -LI d (0 ) = B (c e -d e ) az II o o -3/2 LI -LI - 1 6 (c e +d e ) d&, (4.10) 2 o o az where zb = / Sdz. LI Ztl By use of the WK3 assumption in Eg. (2.27) , drop the second term in Eg. (4. 10), and obtain LI -LI LI -LI y = 6 (c e -de ) / (c e +d e ), (4-11) f o o o o where 6 is evaluated at the bottom of the water column. Solve for c and d , and substitute into Egs. (4.7) and o o (2.26) : LI -LI tan 6 = b /a = £2e 3- (4.12) Note that as LI becomes larger, 9 approaches tt/4 radians or 45°, which agrees with the phase shift for a purely refracted wave in a boundless medium. As LI approaches 64 SOUND SPEED FIGURE 12 - Bottoi Boundary Condition - Case I 65 zero, the turning point approaches the bottom and Q approaches 71.57°. Differentiate 9 with respect to k to obtain 1 r 2 2-1 d 6 = (a +b ) [a b» - b a«], (4.13) dlCpl o o o o o o where primes denote partial differentiation with respect to k . For this case a* and b* are evaluated by r o o LI -LI LI -LI a1 = b Ll» ♦ (2e +e )S ■ - (2e -e )y ■ , o o f LI -LI LI -LI b» = a LI* ♦ (2e -e )g • - (2e +e )y«. (4.14) o o f Ll* and g1 are evaluated by H -1 = d ( edz) = k 6 dz, (4.15) ale J r LI ale ztl and 2 2 1/2 -1 8* = d (k -k ) = k 6 . (4.16) f aTcr r r f The bottom model used must supply a value for u*. 2. Case II For k sufficiently close to k , the wave number at r f the bottom of the water column, the Airy solution (Eg. (3.5)) is a good approximation (Fig 13). From Eg. (3.5) 66 SOUND SPEED :\a,, WKB \psci/latory \Solution LI* Airy Solution FIGURE 13 - Bottom Boundary Condition - Case II 67 2/3 y = [c (/IB -A ) ] / [c (A +/3B )y + 2 o o 1 o o c (/3B -A ) 6 ]. (4.17) 2 oof Solving for A and B gives o o /-, 2 ,_ 2/3 b /a = [c (y^-l) (Y-0 )-c (/3Vl)y Y ]/ O o 2 f 1 2 2/3 [C2(/T+1) (Y-S J-C^/^-IJy Y ]. (4.18) For a rigid bottom y= 0, at k equal to k , 9 equals it/12 r f 1 radians or 15°. If approaches infinity (which corresponds to a pressure release surface), Eg. (4.13) becomes Eg. (4.6) . For is evaluated by Eq. (4.6). For case II the derivative of 9 (Eq. (4.13)) 1 ^ 2 _. 2/3 _ a» = -2c k (/3%1)y - [c 8 (/I+1)+c V G/^-DJv'r o 2 r 2 f 1 b« = -2c k (/3-1) u - [c 6 (/3-1)+c y (/3+1) J y» . (4.19) o 2 r 2 f 1 3. Case III The simplest case is when the wave is reflected from the bottom at grazing angles great enough to allow use of the WKB oscillatory solution (Fig 14). The derivative of the WKB oscillatory solution (Eg. (2.24)) is k dz-3 ) , (4.20) z 1 z where the terms for the derivative of k have been dropped z 68 SOUND SPEED bottom Ill O JJ,, WKB ^Oscillatory {Solution I dil , Udl**1 FIGURE 14 - Bottom Boundary Condition - Case III 69 as in Eq. (4.10) » Thus y = k tan ( Z Z I k dz-$ ) ; J z 1 evaluated at the bottom, the expression for 9 becomes 9 = arctan (-u/k ) • (4-21) 1 z j bottom The derivative with respect to k is r 2 2-1 y«) = (y +k ) [-k y»-yk/k]. (4.22) cflcr 1 z z r z ** • The General Cases As with the surface boundary condition there are three ranges of phase speed for which each of the forms of the solution is a valid representation (fig 15) . At relatively low phase speeds, the shadow zone HKB solution is a good approximation (case I) . At phase speeds near the sound speed at the bottom, the Airy solution is a good representation (case II) . At high phase speeds near cut-off, the oscillatory WKB solution is a good representation (case III) . In the case of the bottom boundary, however, there are two gaps between the representative solutions: one between case I and case II, and one between case II and case III. In a method similar to the surface case, the gaps between these cases can be bridged with two third-order polynomials in k . r C. BOUNDARY CONDITION TEST The equations for the bottom boundary condition were tested for three test cases by comparing the WKB results (as 70 80- 60- 40- 20- CASE I WKB Exponential Sound Speed at the Bottom I "GAP" CASE gap' CASE. WKB oscillatory^ PHASE SPEED FIGURE 15 - Bottom Boundary Condition -9 vs c 1 P 71 calculated by the program QMODE) with a finite difference technique. The finite difference technique was adapted from a normal mode program by the author [ 7 ]. The routine started with a fluid bottom boundary condition and iterated a solution for the depth function upward over a depth grid of 501 points using the finite difference equations of the N0BM01 program of Ref. [7]. At the grid points the sound -2 speed profile was interpolated to be linear in c . This was done to insure the closest possible match between profiles. The iteration was stopped when the lower turning point had been crossed and -2 jk (i+1)-k (i) | h k < 0.1, (4.23) z z z where h is the grid step size. This condition allows the WKB form of the solution to be a fairly good approximation. The lower turning point phase angle was then calculated using the expression -J Q = arctan [-(Jf/(* 0) ] - I * dz 1 z J z ztl where U1 is the derivative of 0 with respect to depth (this quantity is a part of the N0RM0 1 iteration scheme). The integral of k was calculated using the trapezoidal rule. z If Eg. (4.23) was not met, no comparison was made. The three test profiles consisted of two deep ocean profiles and one shallow water profile. The three profiles were the three test cases for the Acoustic Environmental Support Detatchment Workshop on Acoustic Propagation Modelling by Non-Ray-Tracing Techniques [8]- The profiles are listed in Tables 1 through 3. The results are plotted begining with Fig 16. 72 D, REFERENCES 1 Brekhovskikh, L.M., Waves in Layered Media, Academic Press, 1960, p. 2 12. 2 Schiff, L.I.., Quantum Mechanics, McGraw-Hill, 1955, p. 190. 3 Tolstoy, I. and C.S. Clay, Ocean Acoustics, McGraw-Hill, 1966, pp. 48-50. 4 Schiff, L.I., ibid., p. 190. 5 Tolstoy, I. and J. May, "A Numerical Solution for the Problem of Long Range Sound Propagation in Continuously Stratified Media, with Applications to the Deep Ocean", J. Acoust. Soc. Am., Vol. 32, p. 655, June 1960. 6 Bartberger, C. and L. Ackler, Normal M.ode Solutions and Computer Programs for Underwater So. unci FropaqaTionT Part T, Haval~ Air Development Center Report Ho. NADC-72001-AE, 4 April 1973. 7 Evans, K.E., Two Methods for the Numerical Calculation °f Acoustic Nor mal~M"ocIes in th"e 0"cean, H.ST" Thesis, Haval Postgraduate ScEool, Sept 1973. 8 Spofford, C.W., A Sy_no2sis of the AESD Workshop on Acoustic Propagation Hodellinq by ""TTon-ETay-Tracinq Techniques, A~5ED" TecEnical Note TN-/3-U5, N"ov TT73 , p.2"2T~ 73 1 1 • 1 • i t OMODP surface • .«■ '"IWKB eu QMODE bottom m NORMO bottom-— 1 mmmm O co - ■ I f i l ■ / CO \ 2 <3> 00 c CM JO Q UJ UJ m 0. CM* CO Q UJ CO o CO < CM X 1 IO 0. 63 O CM IO 00 O IO 'e puone 'S33a93a "! 319NV 3SVHd 75 IO o CO ro 10 eg 09 ro > IO ex ■"• o m a. o CO t m V. ^ r 00 <3> CM c IO O (^\ UJ CO W LL << CO 1 Tf id CO CO .— 1 CVJ < IO X 0. w o O CM IO 001 08 09 0* r 02 319NV 3SVHd CO IO 76 TABLE I AESD 1 PROFILE AND RON TIMES PROFILE Depth meters 0.0 152.40 406.30 1015.90 5587. 89 Sound Speed m/sec 1536.50 1539.24 1501. 14 1471.88 1549.60 Depth feet 0.0 500.00 1333.00 3333.00 18333.00 Bottom Sound Speed: 1555.52 m/sec. Bottom Density: 1.9176 gm/cc Source Depth: 253.90 m. Receiver Depth: 863.50 m, Frequency: 25.00 Hz Sound Speed ft/sec 5041.00 5050.00 4925.00 4829.00 5084.00 5103.42 ft/sec 833.00 ft 2833.00 ft QMODE RUN TIMES Calculation of 44 Modes: 3.62 sec Transmission Loss to 120.0 nm 18.43 sec 77 TABLE II AESD 2 PROFILE AND RUN TIMES PROFILE Depth meters 0.0 60.96 91 .<*4 182.88 487.68 1219.20 1402.08 1828.80 2438.40 3962.40 5486.40 Sound Speed m/sec 1519. 12 1520.34 151 1.20 1507.24 1503.88 1511.81 1508.76 1499.62 1504.80 1527.66 1554.48 Depth feet 0.0 200.00 300.00 600.00 1600.00 4000.00 4600.00 6000.00 8000.00 13000.00 18000.00 Bottom Sound Speed: 1560.39 m/sec, Bottom Density: 1.9176 gm/cc Source Depth: 243.84 m , Receiver Depth: 1097.23 m, Frequency: 50.00 Hz Sound Speed ft/sec 4984.00 4988.00 4958.00 4945.00 4934.00 4960.00 4950.00 4920.00 4937.00 5012.00 5100.00 5119.40 ft/sec 800.00 ft 3600.00 ft QMODE RUN TIMES Calculation of 81 Modes: 4.36 sec Transmission Loss to 300.0 nm 21.60 sec 73 TABLE III AESD 3 PROFILE AND RUN TIMES PROFILE Depth Sound Speed Depth Sound Speed meters m/sec feet ft/sec 0.0 1524.00 0.0 5000.00 15.24 1520.95 50.00 4990.00 45.72 1517.90 150.00 4980.00 Bottom Sound Speed: 1541.31 m/sec, Bottom Density: 1.40 gm/cc 5056.80 ft/sec Source Depth: 6.10 m, Receiver Depth: 12.19 m. Frequency: 500.00 Hz 20.00 ft 40.00 ft QMODE RUN TIMES Calculation of 5 Modes: 1.32 sec Transmission Loss to 25.0 nm 6.05 sec 79 V. THEORY: THE INgOMOGEMEOOS EQUATION In the preceding section approximate solutions to the homogeneous differential equation (2.10) were developed. In this section a source term will be added to the differential equation, resulting in an inhomogeneous differential equation to be solved by a Green's function superposition of homogeneous solutions. A. SOLUTION OP THE INHOHOGENEOUS EQUATION The inhomogeneous version of equation (2.10) for a source at depth z , and zero range, is s J_ d_(rd$) + d2 $ + w2 = -2 5j[rl 5 (z-z ), (5.1) r "3r "3r Hz7 c? r s where 6 is the delta function. Applying the Hankel transform (Eg. (2.11)) results in the inhomogeneous differential equation 2 2 d*a + [k - k ] u(z,k ) = -2 5 (z-z ). (5.2) ctz? r r s The conditions imposed upon this equation are the same as (2.14) with the exception that du is discontinuous at depth d"z z=z . First, the Green's function will be defined and then s the solution will be integrated using the inverse Hankel transform,. 1 . The Green1 s Function Equation (5.2) can be solved by means of a Green's 80 function £1], defined by u = G(z,z ) = -20 (z) a (z ) / W(0 ,U ,z ) s 12s 12s 0 < z < z , (5.3) s and u - G(z,z ) = -20 (z) 0 (z ) / W{0 ,0 ,z ) s 2 1s 12s z < z < z . s b The function 0 (z) is the solution to the homogeneous equation (2.13) which satisfies the surface boundary condition. The function 0 (z) is the homogeneous solution 2 which satisfies the lower boundary condition. W (0 ,0 fz) is 1 2 the Wronskian for the two solutions, given by W(0 ,0 ,z) =0 (z) 0 • (z) - 0 • (z) 0 (z) , 12 12 1 2 where primes denote partial differentiation with respect to depth. The separate solutions 0 (z) and 0 (z) become 1 2 -V2 f 0 (z) = k cos( | k dz - 0 ) , (5.4) 1 z „J z 1 ztu -1/2 ptl 0 (z) = k cos( J k dz - 0 ) , 2 z J z 2 where z < z < z ; tu tl and 0 and 0 are determined by the surface and bottom 1 2 boundary conditions. Designate ■,-/ k dz - 9 , (5.5) 1 tu 81 /ztl L = J k dz - 9 , 2 r-k (k ) dz - d (« +« ) ]. (5.8) il 2 r J z aiJ 1 2 Ztu r Since W1 is well defined at the zeros of the Wronskian, the poles of the Green's function are simple poles. 2- Evaluation of the Integral To evaluate the integral expression for (z,r) , r $ (z,r) = I G(z,z ;k )J (k r) k dk , (5.9) 0J srOr rr it is standard practice to express J (k r) as a pair of o r Hankel functions (Ref. [2], [3], [4]) 82 ( 1 > < 2) J (k r) = 1 [ H . (k r) + H (k r) ]. or 2 ° r ° r Now becomes (1) (z,r) = 1 I G(z,z ;k ) H (k r) k dk 2 OJ srO r r r J G(z,z ; r (2) 1 G(z,z ;k ) H (k r) k dk . (5.10) 2 oJ srO r rr + 2 Two properties of the Hankel functions are now used: (1) H ■+■ 0 for large argument in the first quadrant of the 0 (2) complex plane, and H -+• 0 for large argument in the fourth quadrant. The first integral is evaluated by integrating around the first quadrant, while the second integral is evaluated around the fourth quadrant (Fig 19). The poles of the integral lie on the real axis so the integrals in (5.10) are not well defined. However, the boundary conditions are met only by integrals whose path of integration in the complex plane passes below the poles. The real part of the horizontal wave number, k , represents the propagation of r the wave functions with range. Any imaginary part of k r represents an expontential growth or decay of the wave function with range. For k to represent an outgoing wave r (taking into consideration the sign in Eq. (2.8)), any negative imaginary part represents waves which "exponentially grow" with range. Since this is inconsistent with the original restrictions on T , it is required that 83 imaginary Complex kr Plane I 03 Branch Point Vi ROLES, krn a'2 kr w/cb FIGURE 19 - Path of the tfankel Integrals in the Complex k ] Plane 84 Im(k ) > 0- This is physically appealing because it is r equivalent to having a small amount of attenuation. In addition, there is a branch point at k equal tou/c , where r b C is the sound speed in the bottom. However, at long range b the associated branch line integral has been found to make an insignificant contribution [5]. The first integral now may be represented by ioo \ a0 fix 0 and the second by ) (1) / H" (k r) G k dk ♦ TriZRes{H (k r) G k }, (5.11) 2^r0r rr Or r o"i J (2) H (k 3 r2 0 r 1 / H (k r) G k dk . (5.12) 2 Jrr 0 r r r -ioo Since (1) (2) H (ik r) = - H (-ik r) Or Or and G(z,z ;k ) is an even function of k , the integrals in s r r equations (5.11) and (5.12) cancel. This leaves 8 (1) (z,r) = +Tri I Res{H (k r) k G(z,z ;k )}. (5.13) 1 0 r r s r Since the poles of G(z,z ;k ) are simple and occur at the s r zeros of the Wronskian, the derivative of the Wronskian with respect to k is used to find the values of the residues. r Further, since at the zeros of the Wronskian, the two solutions U (z) and U (z) are identical let 85 U = (J = 0, 1 2 which yields (1) (z,r) - -2iri S [H (k r) U(z)0(z )k ] / W« . (5.14) poles Or s r 3 . Normalization From the derivative of the Sfronskian (Eg. (5.8)) a normalization factor is designated ftl I -1 -1 ||u2jj = I k dz ♦ k d f fl + e ]. (5.15) J z r elk" 1 2 Ztu For long range, equation (5.14) can be written (using the asymptotic form of the Hankel function [6] 1/2 -1/2 i Q = 1 [ kdz-e-6+irl, 11 Z ,J z 2 1 ir z~ and the minimum Q of the main channel is ■ i '/" Q = 1 [ I kdz'-e - « +tt]. Ill ■„ ZJ Z U 1 tu The above expressions for Q , Q and Q show that the sum I II III of maximum modes in the upper and lower channels may be equal to N-2, N-1 or N. Thus it is possible for the WKB characteristic equation to have no solution or multiple solutions for a particular mode near the crossover between multiple channel and main channel families. This phenomenon is a mathematical consequence of the fact that near the edge of the barrier, both k and g are sufficiently small to z invalidate Eq- (3-17). This difficulty could be overcome in a number of ways, which "include changing the nature of the barrier or using a method other than the WKB method in the region near the barrier edge. If the edge of the barrier is appropriately smoothed so that the derivative of the vertical wave number is proportional to the square of the vertical wave number, the WKB validity criterion will remain satisfied. For this particular model this scheme would be undesirable, as it would require the use of a profile other -2 -2 than the c linear profile. The c profile is used primarily for the rapidity with which the vertical wave number may be integrated. Since integration of the vertical wave number is central to many operations in the WKB method, this integration time is extremely important to the program execution time. Another method to overcome this difficulty 98 would be the alteration of the barrier connection formulae (Eg. (3.15)) for very small values of the integral L. A third method would be to modify the integral of the vertical wave number, so that a correction term is applied at depths which fail to meet the WKB validity criteria. Again, this method adds to the execution time required for the integration loop. In the program described in the next chapter the difficulty of a duplicate eigenvalue is overcome by eliminating the duplicate eigenvalue from the main channel. The difficulty of a missing eigenvalue is overcome by assuming the missing eigenvalue to be at the horizontal wave number corresponding to the edge of the barrier. B. EIGENVALUE ESTIMATION Following definition of the mode families, the next step is the determination of the mode eigenvalues. 1 . The Eigenvalue Search The search for eigenvalues is done by families, determining those values of the horizontal wave number which most nearly yield integer values of Q. The eigenvalues are found in three steps; the formation of an initial polynomial estimator; the finding of two values of k which r bracket k ; then performing successive iterations by the r,n method of regula falsi. For each family, sample values of Q as a function of phase speed are found (by integration of the vertical wave number and determination of d and $ ) . 1 u Through these points a fifth-order polynomial in Q(k ) , P (Q) i is fitted by the method of least squares. This 99 polynomial is used to find the first estimate of the horizontal wave number. For long ranges and when considering incoherent transmission loss, this estimate is sufficiently accurate, and the subsequent iterations may be bypassed. Using the derivative of the polynomial estimator, p» (n) , a second estimate for k is then sought which, in r,n conjunction with the first, brackets the eigenvalue. From these two estimates which bracket the eigenvalue, the eigenvalue estimate is improved by successive iterations, up to four times (the method of regula falsi [ 1 ]) . 2- Selection of the Estimator A number of techniques were considered as estimators for k as a function of Q, including a cubic spline and the r least squares polynomial (which was used) . a- Cubic Spline The cubic spline fits a piecewise cubic polynomial through a given set of data, resulting in a function continuous through the second derivative. When Q(k ) was monotonic and reasonably smooth the spline gave r excellent results (in fact, it was exact for certain analytic profiles) . In the case of the multiple channel profiles, Q (k ) may become irregular and slightly r oscillatory in character. In this case an unfortunate combination of sample points will generate a spline which gives estimates of k in the wrong order (Fig 24) . The r cubic spline was deemed too sensitive to the minor phase effects of the multiple channel modes, and a smoother substitute was sought. 100 Y vS,. \ \ » r ! 1 1 * " * • «^*_ X, » Xx * ^. % ^>^ 0 sample points — cubic spline • r— oft,; II Q FIGURE 24 - Ill-conditioned Cubic Spline 101 brf Least Squares Polynomial The problems with the cubic spline led to an estimator which would be insensitive to any small scale oscillation in the data, yet which reliably estimated the general trend of Q (k ) . A least squares polynomial r algorithm using orthogonal polynomials from the text by Conte and deBoor [2] was selected. To decrease the difficulties near the edges of multiple channel barriers, a weighting was used which de-emphasized the minimum and maximum values of Q in the families. The polynomial was found to be of sufficient accuracy to provide, in most cases, an estimate accurate to within one part in twenty of the intermode spacing. In addition, the polynomial provides a rapid estimate of the derivative of Q(k ) with respect to r k ; this result is useful in later iteration steps and in r normalizing the eigenf unctions. C- CALCULATION OF THE EIGENFUNCTION After the eigenvalue ic is determined, the next step r,n is to calculate and normalize the eigenf unction U (z) at the n source and receiver depths, z and z . To determine the s r eigenf unction, one of the three solutions developed in the theory chapters {the WKB and Airy solutions) , or a function extrapolated from the known solutions is used. 1 . Assumed Solutions At the receiver and source depths, the Airy and WKB validity criteria, F and F , are evaluated. If a WKB or w a 102 Airy solution is a good approximation, the eigenf unction is determined in two steps. ■ The coefficients for the WKB solutions, a , b , c , o o o and d are determined for each ensonified and shadow zone o region, using the values of the upper and lower turning point phase angles d and 0 . u 1 ■ If the WKB solution is a good approximation, the integral of the vertical wave number from the upper or lower turning point to the source or receiver depth is determined, and the eigenf unction is evaluated using Eg. (2.22) or (2.26). ■ If the Airy solution is a good approximation, the Airy coefficients, A and B , are determined from the o o connection formulae (Eg. (3.13)); the variable x is calculated; and the Airy solution is evaluated with Eg. (3.5) . 2. Spliced Solution If neither the Airy or WKB solution is a good approximation at the desired depth, a finite difference technique is used to extrapolate the eigenf unction from the two accurate WKB or Airy solutions above and below the desired depth. A search is made from the desired depth to find the nearest depths above and below, z and z , for a b which either the Airy or WKB forms are good approximations (Fig 25). At these depths U (z ) and 0 (z ) and their a b derivatives with resDect to depth, U * and U ■ , are a b calculated. The depth interval between z and z is divided a b into (up to 100) equal increments, z , and a solution is i 103 3* t T v. |N -p--' Q. UJ co "to la to tu Z. 5 o Q. N N A. N ^ U. N 2 tt / X . r,nrm s A running average, which is a window in range space, represents a filter in wave number space. Thus a running average of the acoustic pressure squared term, optionally weighted with a Gaussian filter, is used to give an averaged transmission loss figure. The Gaussian filter was selected since the Gaussian weighting function in ran-ge transforms to a Gaussian weighting function in wave number space, allowing a smooth roll-off from long wavelength to short wavelength fluctuations. This physically corresponds to the elimination of the finer details of the sound speed profile as contributors to the calculated acoustic field. 3. Incoherent Transmission Loss At long range ("long" depends upon the acoustic frequency and the inhomogeneities in the ocean) , the coherent interference effects lose significance and the incoherent acoustic pressure (representing an SMS sum of the energy in the modes) becomes the most useful estimate for the acoustic field. Accordingly, incoherent transmission loss is defined 2 TLI = -10 log Z P . (6.10) 10 mm E,. MULTIPLE PHOFILE PROPAGATION In the preceding sections no horizontal variation in both sound speed and bottom depth has been assumed. Within one locale (at sub-kilohertz frequencies) this assumption 108 can be fairly accurate. However, as sound propagation over greater ranges is considered, markedly differing oceanic regimes are encountered, and this assumption becomes untenable. Thus the long range sound propagation model must include the effects of changes in the sound speed structure with range. 1 • Adiabatic Assumption To incorporate horizontal change the so-called adiabatic assumption is invoked, namely: the energy in any particular mode stays within that mode. The assumptions involved may be more explicitly stated. To paraphrase Tolstoy and Clay [4]: ■ The eigenf unctions and eigenvalues are dependent upon the local stratification. ■ The acoustic medium varies slowly with range. ■ There is no cross coupling between different modes in the transition between profiles. ■ There is no back reflection of acoustic energy at the transition between profiles. 2- Adiabatic Transition Range Milder has shown that the accuracy of the adiabatic assumption decreases with increasing frequency and 2 (a) (b) horizontal gradients of k [5]. If u and u are m n adjacent modes (nm±1) , Milder shows the adiabatic assumption is a good approximation if JZ;1 (a) 2 (b) 2 2 u d (k ) u dz} / k « 4tt , (6.11) m aT n where X is the skip distance for an equivalent ray. Milder1 s equations may be extended to estimate an order of magnitude transition range required for an adiabatic change 109 between two profiles. If the magnitude of the integral in Eg. (6.11) is estimated by z r1 (a) a ,v2. (b), I u d (k ) u dz J m -» or r 2 • 2 ^lldl*-)-! dz = yj -i- ,Ak ' dz' (6.12) and an expression for Ar, the distance required between the two profiles for an adiabatic change to occur, is found 2 2 2 Ar » X IaJc Idz / (4/2 tt ). (6.13) This transition range is calculated by the program and printed out as a guide to the applicability of the adiabatic assumption. 3- Transmission Loss Since the eigenf unctions and eigenvalues are assumed to be characteristic of the local profile, the source and receiver mode excitation is expressed in terms of the (s) (r) eigenfunctions at the source and receiver, u and u n n respectively. Defining : (r) = / k dr, (6.14) n J r,n and D (r) =| D dr, (6.15) n J a 0 the adiabatic mode pressure term is derived, 1/2 (s) (r) P = 2[2 tt/K (r) ] u u exp[-D (r) ]. (6.16) m m mm m This gives the transmission loss for the multiple profile case: 110 2 TL = -10 log [ Z P cos(K (r) ) ] + 10 n n n 2 [ IP sin(K (r) ) ] . (6.17) mm m 4. Mode Matching Multiple profile propagation using the adiabatic assumption requires the matching of modes between profiles. In previous works this has been done by numbering the modes in order of increasing phase speed. This is always correct for the single channel profile, in which the mode order is always invariant. However, with multiple channel profiles this procedure may yield physically unreasonable results. Consider two sound speed profiles with two near symmetric channels, and six modes (Fig 26) . Since the upper channel is slightly "deeper", the modes in this channel occur at a lower phase speed than the corresponding modes in the lower channel. Thus, in terms of overall order, the odd-numbered modes represent acoustic energy trapped in the upper channel and the even numbered modes represent acoustic energy trapped in the lower channel. Slowly alter the profile so that the lower channel becomes slightly "deeper" than the upper; and the modes in the lower channel precede the modes in the upper channel. Matching modes by overall order would, in this case, couple all the modes in the upper channel with modes in the lower channel, and vice versa - a physically implausible event. In matching between profiles, not only mode order, but mode families must also be taken into account. Assume again the profile changes slowly, this time so that both upper and lower channels become "shallower" (Fig 27) . Two modes appear in each of the upper and lower sub-channels, and two modes appear in the main channel. It is not clear which of the new modes in the main channel corresponds to the missing last modes from the upper channel 111 V) O a -p •H « 0) c c (d U o •H M -P 0) a 3 >i M (U 25 o 3 EH CM W M O H Eh Wop^ 112 1 + F (Q/2* until r a value for k<2> is found for which, r,n (Q(k«))-n)*(Q(lt(l))-n) <0. r,n r,n Once two bracketing values are found, the iteration continues using the method of regula falsi until one of two conditions is met: a value of Q within 0.0005 of the integer n is found; or a total of five iterations have been made. Once the eigenvalue k has been found, a number of r,n mode parameters are then calculated and stored for later 122 use. These include: k (KRVEC) , mode order (MODENO, r,n MORDR2) , the mode number (MODEN) , relative position in the sound profile (IBRN2) , the skip distance (SKPDIS) , and Q (k ) (QN) . The subroutine EIGFDN calculates the value of r,n the eigenfunctions at the source and receiver depths. The eigenf unctions are then normalized, completing the mode and family loops. If the current profile is the continuation of a multiple profile run, the current mode set is matched with the source location modes by the subroutine MATCH. The subroutine PLOSS calculates the transmission loss from the range RSTART to RSTOP in increments of DELR. If the next profile is to be a continuation of the same run, two subroutines are called to prepare for the continuation. The subroutine PLSET attenuates the source eigenfunctions and keeps a running sum of K (r) (KRRV) ; subroutine SWAP keeps track of n the matching parameters. If, on the other hand, the next profile is the beginning of a new run, the default parameters are re-set. The input data are then read for the next profile and the main program sequence begins again. B. SUBPROGRAMS 1 . BOTSET The subroutine BOTSET selects the bottom type (IBT) , sets parameters for the calculation of the bottom boundary parameter, y , and the mode attenuation. The functions EMU and EMUP compute y and y1 as a function of k . The function DECAY computes the mode attenuation as a function of the bottom type and characteristics, bottom roughness 123 (BOTfiUF) and surface roughness (WAVEHT) . If the modified fluid or the modified rigid bottom models are selected, DECAY linearly interpolates the bottom loss values from input data of bottom loss vs grazing angle. 2. BCSET. The subroutine BCSET divides the bottom and surface boundary conditions into phase speeds based upon the values of the Airy and HKB validity factors, F and F , versus the a v parameters KALIM and KWLIM. For the surface boundary the HKB forms are used at phase speeds below CP3 and above CP1 ; and the Airy forms are used at phase speeds between CP1 and CP2. For the bottom boundary the WKB forms are used at phase speeds below CP7 and above CP4; the Airy forms are used at phase speeds between CP5 and CP6. For the "gaps" between these phase speeds (two at the bottom boundary, one at the surface) , third-order polynomials are formed to provide continuous functions for -9 and 9 . The subroutines u 1 SBC and BBC calculate the values of 9 and 9 as functions u 1 of 3c . In order to distinguish between multi-channel r families, a minimum or maximum family depth is required for SBC and BBC respectively. When the parameter JP is set to -1 the subroutines also calculate values for 9 ■ and 9 • . u 1 3. FIN F AM The subroutine FINFAM defines each mode family by checking for family boundaries at the surface and bottom sound speeds, and each relative sound speed maximum. Two arrays are used to store the parameters for each family: maximum k , minimum k , minimum Q, and maximum Q (in FAM) ; r r and minimum depth grid, maximum depth grid, family type and 124 relative position in the sound speed profile (in IFAM) . 4. ORTPOL This subroutine, from the text by Conte and de Boor, [2] forms a fifth order least squares polynomial through a set of data points, representing values of Q (k ) and k , r r each with a weight W. The resulting polynomial, evaluated by the subroutine ORTGET, finds the first estimate of Jc r,n for a given value of Q. ORTGET returns a value of the polynomial and its derivative in terms of both phase speed (CPV) and horizontal wave number (KM) , as well as the derivative of the k with respect to Q (DKRDQ). r 5. EIGFON Given a value of k and a depth, z, (entered by r,n depth grid index (IDZR) and the depth below the grid depth (DELZR) ) , the subroutine EIGFUN calculates the unnormalized eigenf unction U(z)» Upon the first call to EIGFON after the determination of the eigenvalue, the WKB coefficients, a or o c (C01) and b or d (C02) are calculated for each depth o o o grid as a function of the upper and lower turning point phase angles , i may exponentially grow with depth. In order to force U to i be continuous at z , u. This in turn is multiplied by a second ramp function i yielding the extrapolation function U<*>, which converges to i U (z ) . There are now two continuous extrapolation functions a U<2> and U<4), with U<2> representing a more accurate i i i extrapolation near z and U <♦) representing a more accurate a i extrapolation near z . An average of the two extrapolators b weighted by the relative proximity of the desired depth to z and z is used for the final estimate of the a b 126 eig en function, U (z) = [ (z-z ) U<2> + (Z -z) U<*>] / (z -z ) . n b i a i a b Due to the ramp function used, this extrapolation scheme may give a discontinuous derivative at the end points, z and a z . Since the derivative of 0* (z) with depth plays no part b in the calculation of the transmission loss, and the errors incurred in 0(z) are slight, this effect can be expected to have little if any effect on the subsequent calculations. 7. PIjOSS The subroutine PLOSS computes the transmission loss with range from RSTAHT to RSTOP in increments of DELE, and displays the results on a printer graph or a CALCOMP plot. The transmission loss is calculated in terms of both the coherent and incoherent transmission loss, and an optional averaged transmission loss. The averaged transmission loss is calculated from a running average of the square of the 2 coherent pressure P multiplied by a range weighting function (or window) . Two windows are provided, a Gaussian curve (IFILT=1) and a straight running average (IFILT=2) . When PLOSS is called, various initial values are set and the filter function is normalized. The averaging is 2 performed with a stack which is filled with P values. At each range step the values in the stack drop one place; the 2 value of P is calculated for the range r plus the averaging width (RAVG) and is placed at the end of the stack. The coherent transmission loss is calculated from the middle value of the stack. 127 For multiple profile plots the transmission loss calculations stop at the range equal to the end of the profile sector (REND) , and a variable is set indicating the subroutine was left "on" ($0N) . After processing of the next profile, the transmission loss calculations continue without interruption to the output. 8. MATCH This subroutine fills two vectors (ICSOSS ,IREX) which provide a cross-reference table for the matching of modes between profiles. The matching is carried out in two steps: ■ Each mode in the new profile is selected in order of increasing phase speed (M0RDR2) . For each mode, the modes in the previous profile are searched (also in order of increasing phase speed) for a mode with the same relative position in the profile (IBRN1,IBRN2) and with the same mode number (MODEL, MODEN) . ■ If, after the first step, any modes remain unmatched a transition between channels has taken place. In this case the second step begins by finding the extent of the gap of unmatched modes. The gap begins with the first unmatched mode in the new profile, and ends with the next matched mode in either the same channel or the main channel. The corresponding gap is also found for the previous profile, and the number of missing modes in each channel is tabulated. The modes are then cross-referenced (ICROSS) in terms of increasing phase speed. As long as unmatched modes remain in more than one sub-channel, the source mode energy is summed (USUM) , later to be divided evenly between the same modes. This step is repeated for all such gaps. Finally, any modes which still remain unmatched are cross-referenced to mode 200, which is used as a dummy, or null mode. 128 9 - INTEGR The subroutine INTEGR integrates either the vertical real wave number or vertical imaginary wave number with depth. With the assumed nature of the profile this integration is accomplished in each layer with the analytic formula, J V2 = 3 3 2 |[lc (z.) - k (z. J ]| / y 3 z i z i+1 i 1 0 . £F I ND The subroutine QFIND computes the value of Q at a particular k for a family between the depth grids MINDEP r and MAXDEP. This is accomplished by integrating the vertical wave number across the ensonified depths (INTEGR) , then calling the routines for the surface and bottom boundary conditions (SBC and BBC) . C. REFERENCES Wilson, Wayne D. , "Speed of Sound in Sea Water as a Function of Temperature, Pressure and Salinity1', J. Acoust. Soc. Am., Vol. 32, No. 6, p. 641, June 1960. Con Ana nte, S. D., and C. deBoor, Elementary Numerical alysis, 2nd ed. , McGraw-Hill, 1972, p.25Sr26F: 129 VIII. RESSiilS AND CONCLUSIONS Results form the QMODE model were compared with analytic solutions, the results from some other models, and with some transmission loss experimental data. A. ANALYTIC TESTS The model was tested with two profiles which yield analytic values for the eigenvalues and eigenf unctions. 1- Isovelocity Channel The first profile is that of an isovelocity sound channel with rigid bottom, for which the eigenvalues are 2 2 2 k = k - [ (2n-1) it/(2H) ] , (8.1) r,n where H is the depth of the channel. The profile had a sound speed of 15G0 m/sec and a depth of 4000 meters. The first 150 eigenvalues were computed for a frequency of 100 Hz by both the QMODE and NORM01 (a normal mode finite difference program by the author [1]) programs. The results are shown in Fig 28, where the error between the analytic and computed eigenvalues, scaled to the intermode spacing, is plotted as function of mode number. A curve is also included which corresponds to six significant figure agreement between the computed eigenvalues and the analytic eigenvalues. For higher modes, as the vertical grid spacing becomes an appreciable fraction of the vertical wavelength, the finite difference technique becomes more inaccurate. 130 2. Constant G radiant Channel -2 The second profile is one with a single c gradient (Fig 29) . For this profile the analytic eigenf unctions are a linear combination of the Airy functions, Ai and Bi. In the shadow zone 3i increases exponentially with depth, for modes with lower turning point well above the bottom, and the coefficient of the Bi term must be zero. Thus the eigenvalues are those values of the horizontal wave number for which the zeros of the Airy function occur at the surface. This occurs when (using the notation of chapter III), Ai(-x) = 0, (8.2) or, 2 2/3 2 k = Y y - * * (8-3) r,n n o where y are the negative zeros of the Airy function and Jc n o is the wave number at the surface. The differences between the analytic eigenvalues and the computed eigenvalues are shown in Fig 30 as a function of mode number. For the lower modes the finite difference technique and the 8KB model have comparable accuracy. Beginning with mode 12 the QMODE (WK3) eigenvalues become accurate to six significant figures, asymptotically approaching the analytic solution. B. AESD HOfiKSHO? TESTS The QMODE results for three profiles used in the AESD Workshop on Non-2ay-Tracing Techniques were examined [2]. The three profiles correspond to a Pacific deep water profile (AESD 1B) , a N.E. Atlantic profile with multiple sound channels (AESD 2) , and a shallow water profile 131 (AESD 3) . The sound speed profiles for each are shown in Fig 31, 32, and 33. 1 . Comparison of Modes From two of these profiles, a number of eigenf unctions Mere calculated by both the QMODE and NORM01 programs and plotted. All five eigenf unctions for the shallow water case are shown in Fig 34. The normalized eigenfunctions generated by both programs are superimposed, with the N0BM01 modes in solid lines, and the QMODE modes in dashed lines. Modes 11 through 20 of the multiple channel profile, AESD 2, are shown in Fig 35 and 36. These modes occurred at phase speeds near the edge of the barrier separating the two channels. 2. Transmission Loss Transmission loss was calculated by QMODE for the AESD profiles. The QMODE results are compared with those for the FACT (modified ray theory) , the NOL (finite difference normal mode) , and the parabolic equation models. A description of the NOL model is found in Ref. £3]. The central processing unit (CPU) times for each QMODE run are given in Tables 1 through 3. The transmission loss curves are given beginning with Fig 37. The QMODE transmission loss plots may show up to three curves: the fully coherent transmission loss, the incoherent transmission loss, and the averaged transmission loss. The QMODE results agree generally, although not necessarily in detail, with those for the finite difference and parabolic equation models. The FACT results are considerably smoothed and do not reflect the sharp interference effects apparent in the normal mode and parabolic equation models. The FACT model incorporates acoustic energy with higher grazing angles at the bottom (greater than critical angle) than normal mode techniques. Thus the QMODE results show greater transmission loss in the shadow zones, especially for 132 profile AESD 1B. If RAVG, the range over which the QMODE transmission loss is averaged, is increased the results are smoothed and begin to resemble the FACT results. This is illustrated in Pig 49 through Fig 52, where the averaging range for the AESD 1B transmission loss is increased from 2.5 nra to 10 not. C. COMPARISONS WITH DATA The transmission loss predicted by QMODE was compared with the results from two sets of acoustic experiments. 1- PARKA IIA Data The first set of comparisons was with a series of experiments in the North Pacific (PARKA IIA) [4]. The experiment measured transmission loss over different paths for various combinations of frequency, source, and receiver depth. The research platform FLIP, equipped with suspended hydrophones, acted as the receiver site. Explosive charges were dropped along radial bearings from the receiver and the resulting transmission loss was recorded. Four runs were compared, two to ranges of 500 nm, one to 1500 nm, and one to 2000 nm. a. Run 1 In this run the transmission loss at 100 Hz between a deep source (500 ft) and a sound channel axis receiver (2563 ft) was measured out to a range of 500 nm. The results are shown in Fig 53. The QMODE results agree well in the convergence zones and at longer ranges, but show greater transmission loss in the shadow zones between the first few convergence zones. This can be attributed to the fact that QMODE was run without including any bottom bounce energy (since detailed bottom information was not readily 133 available) . b. Run 2 This example measured transmission loss at 100 Hz between a shallow sourca (60 ft) and a shallow (300 ft) receiver. The data were collected in two runs, one outbound from the receiver for a distance of 500 no, and then a run back towards the receiver. The results for the outbound and inbound runs are shown ia Pig 54 and 55. Again the QMODE results agree with the data at the longer ranges and in the convergence zones. The effect of omitting bottom bounce energy is still apparent. c. Run 3 The third PARKA example measured the transmission loss at 100 Hz between a deep source (600 ft) and a sound channel axis receiver out to a distance of 1500 nm. The environmental input for the QMODE computer run was based upon the Fleet Numerical Weather Central (FNWC) climatology for the appropriate month. This climatology is formed by averaging the temperature and salinity fields for historical data. The computer run used a profile based on climatology for every 100 nm along the tracic. The results are shown in Fig 56. For this northerly track the sound channel slowly rises with range and the shallow source becomes increasingly coupled to the axis receiver. d. Run 4 The last example from the PARKA experiments was along an east-west track for a deep (300 ft) source and axis receiver out to a distance of 2000 nm. Along this track the sound channel axis depth remained relatively constant. The environmental data for the computer run was based upon the FNWC climatology with the temperature values corrected at shallow depths for observed AXBT (Airborne Expendable Bathythermograph) data collected at the time of the 134 experiment [5]. The results are shown in Fig 57. For the first 1000 miles the QMODE results agree well with the observed results. After this range the QMODE model fails to predict the observed large variations in the transmission loss. The FNBC propagation loss calculations for the same track using tha RP-70 model (ray-trace) also failed to show the observed variations ['4]. It is assumed that the wide variation of transmission loss is due to differences between the assumed and actual bathymetry along the track. 2- Antigua to Grand Banks The second acoustic experiment for which comparisons were made measured tha transmission loss along a great circle track from a location near Antigua, West Indies to the Grand Banks. The results of this experiment were described in an article by Guthrie, Fitzgerald, Nutile, and Shaffer [6]. The author was present at the receiver site in connection with military duties. Two acoustic CW projectors at frequencies of 13.89 and 111.1 Hz were towed at depths of 104 and 21 meters. The received signal was recorded through a 1100-meter deep, bottom-moored hydrophone off Antigua. The QMODE results show the effects of variable oceanographic conditions along the track. For each source frequency two sets of calculations were made, one with a bilinear profile (Fig 58) , and a set of profiles based upon FNWC climatology at 100-nm intervals. For the low frequency source, Fig 59, 60, and 61 show the results for the bilinear profile and the climatology profiles, as well as the observed transmission loss. The QMODE results show well-defined convergence zones along the entire length of the track, while the observed data show the zones to lose definition after about 1800 km. Figures 62, 63 and 64 show the corresponding results for the high-frequency source. Again, QMODE predicted well-defined convergence zones to greater ranges than was observed. The appearance of well defined modes at greater 135 ranges suggests that more scattering of acoustic energy among the modes took placre than was accounted for in the computer calculations. This scattering can be due to either inhomgeneities in the ocean or the non-adiabatic coupling among modes of different profiles. Guthrie, et al, observed that as the track extends to the north the sound channel becomes shallower and more acoustic energy is coupled from the shallow sources into the deep sound channel. This change in the nature of the sound speed profile is observed in the QMODE results, where the incoherent transmission loss decreases as the northern latitudes are reached. D. EFFECT OF THE BOTTOM ON DEEP SOUND CHANNEL PROPAGATION In the PABKA run 4, lack of information on the bottom appeared to have an adverse effect on the computer results. Due to a lack of data, the computer runs were made considering only refracted energy (this approximation is often used when dealing with such long range transmission) . Along tracks where the modes have little interaction with the bottom, this limitation has little effect. This usually occurs in the Northern Hemisphere for northerly tracks where, typically, the sound channel becomes shallower with increasing latitude. On an east-west track (as in PARKA run 4) , and on southerly tracks (where the sound channel becomes deeper with lower latitudes) , the modes tend to interact more with the bottom. For such tracks the calculation of deep sound channel transmission loss becomes difficult without detailed knowledge of the bottom characteristics. Thus to predict low-frequency sound propagation reliably over long distances in the ocean, detailed bottom information is required in a form suitable for the mathematical model being used. 136 E,. CONCLUSIONS QMODE possesses the potential of a relatively fast method for the calculation of long range low-frequency sound transmission in the oceans. The computational speed results from the normal mode method of computing transmission loss and the method in which the eigenvalues and eigenfunctions (modes) are determined. The theoretical basis for the eigenvalues and the eigenfunctions is the assumed WKB solution. The first estimate for the eigenvalue is obtained by using a polynomial estimator through data generated by the WKB characteristic equation. When a precise eigenvalue is required, this estimator allows the use of significantly fewer iterations. For the calculation of averaged or incoherent transmission loss, the estimator alone provides an accurate enough approximation of the eigenvalue. Since the eigenf unction is an assumed WKB function, a very precise eigevalue is not required in order to generate a good eigenf unction. This contrasts with conditions in the finite difference technique, where a precise eigenvalue is required due to instabilities in the finite difference equation in the shadow zones. The second characteristic of the model which allows rapid computational speed is a property common to all normal mode models: the accuracy of the transmission loss calculations does not depend upon the size of the range step between calculations. To calculate the transmission loss at some range, the following is required: ■ The modes be properly defined at the receiver; ■ The coupling between modes be evaluated in a range dependent environment, ■ The mode attenuatioa be evaluated; • The modes be defined at the receiver. Thus, the steps required in this model are dependent 137 upon the number of intervening profiles and not the range between source and receiver. In contrast, ray-tracing methods must calculate ray propagation over small grid steps along the entire path. k limitation of the mode method is that the number of modes, hence computation time, is proportional to the frequency considered. There are a number of applications for which the model is ill-suited. The WKB method should not be expected to yield accurate results in situations in which only a few modes are present, such as in a sound channel near the cut-off frequency. As has been seen, the few modes near the edge of a barrier which separates multiple sound channels should not be expected to be exact. In the computation of transmission loss, this should not be a limitation if more than a few modes are present. A weakness of the normal mode technique is that rapid variations of sound speed and bottom depth with range are not easily accommodated. A high relief, sharp bottom becomes extremely difficult to model; this is true as well for rapid changes in the sound speed profile. A number of improvements concerning both the theory and the actual program could be made to the model as presently configured. Some method of estimating the effects of a steep-sloped bottom upon the distribution of energy in the modes would improve the model's ability to predict transmission loss o?ar variable bathymetry. Such an estimate would physically correspond to the non-adiabatic coupling between modes; or in terms of rays, the ray deflection after striking a sloping bottom. Secondly, an estimate of the scattering of energy in the modes due to both inhomogeneities in the environment and to the non-adiabatic transition between profiles would be desirable. A number of programming improvements could be made, including more efficient coding, space for storing a larger number of modes (presently 200) , and multiple source/receiver depths. 138 F. REFERENCES Evans, K.E., Two Methods for the Numerical Calculation of Acoustic ""ITormal "ili^es in tEe""Ocean, M75". TTie*sis, UavalTos^graduale School, Sept T973. Spof ford# C- W. , A Synopsis of the AESD Workshop on Acoustic Propagation' ffoaelTing" 5y TTon-ffa y-*Tr"acin5 Te^Fni g uesTT SED~T e cE n i caI""lTof I TN-73-057~TTo7tT73: Blatstein, Ira M.r Comparisons of Normal Mode Theory, Ray Theory, and Modified Hay Theory"! or~Ar5rErar7~"?o una Velocity Profiles Hasulting in Convergence~Zones, Naval OrcTnance UaForaTor?^"'TecHnicai 'Report 7*4 r95 , 29 Aug 1974. Maury Center for Ocean Science Report 004. PARKA II, A B£i§fi!i2 Report (U) , November 1970, (CONFIDENTTATf. ~ Maury Center for Ocean Science Report 6, Vol. 2, PARKA Oceanographic Data Compendium - PARKA IIA, May 19777 Guthrie, A. N.# R. M. Fitzgerald, D. A. Nutile, and J. D. Shaffer, "Long-Range Low-Frequency CM Propagation in the Deep Ocean: Antigua-Newfoundland", J. Acoust. Soc. Am., Vol. 56, No. 1, p. 58, July 1974 139 T c UJ o t0 a o 2 or. 1 o i 2 o O z jt \ \ I I H 0) a s -C o e > a o> •H w as H Uj pa jotuj^3 FIGOEE 28 - Eigenvalue Test: Isovelocity Channel 140 ROD 1000 M/SL.C TO SOUND SPtE'D URLUE5 S'yjW !P£FJ] JN METFJ^5EC 10. ac J> C3 C) u r- "0 — i T *-4 -.i ir! 71 ^1 .a Lft . o <.) o O FIGURE 29 - Constant Gradient Channel 141 t-U*J U*J, U'J>| 3U/|oud-u% patDUiiisal H 0) a c o o cr LT\ O -P en OJ 3 CO Qi a: O < o 2: 1 a: 0 W « D C5 H fe iG.oc 2J.CC 30. CC DCPTH IN METERS - 1 /-"■ 'Ju • IK i.KIO M.OC 147 IO.vDC CCPTH IN "F.TF.R5 «K10 I 40. CG 50. QC (/> +-> 3 V) O) cc o C\J LU O sz o CT> SL 3 cr O s- V) to O I 10.00 to v> o> -o o 2: CM a LU <: 1 ^o CO LU C3 20. CO 30. CG CF.PTH IN ME1FRS 40. 0C 50. ec 'W0 1 143 QP 'SSOH NOISSlWSNVdl 149 CO 4-> CO cc Q CO LU CO a: qp 'SS01 NOISSIWSNVdl 150 Q CO UJ CO LU cm r> qp 'SS01 NOISSIWSNVdl 151 o 00 .o Ul to CD o en (0 c 0 •1- +-> eo 3 cr O to 0 •r- O £ c o «» 0. 0 OH 0 to r\J 1— 1 .0 qp *SS01 NOISSIWSNVdl 152 +-> HI Q O cr o LU ^~ IO 3 CM CO ~ ~ a) C£ E c C_> ■h < uj Ll_ o& oz 1 — < CM Q CO UJ < lO 1 h- CM «3- LU OC Z2 O t— * O u_ IO q? 'sscn NOissiwsNvyi 154 to •4-> in O) I CM O LlJ CO «3- W 4SS01 NOISSIWSNVdl 155 o o to cc c o •r- =3 IO cr C4 LU u E •r- c o ** -Q LU (0 o o U- o < Q_ 1 CO Q to m LU N- 1 LU CtL o =D IO C3 qp *SS0"1 NOISSI^SNVdl 156 +-> t/j O) Q O cr i CO Q in gp 'ssoi noissiwsnvhi 157 o -«o VI +J 6 3 c C3 gp 'ssoi NOissmsNvyi 159 4- 7.GC d.G : K 10 ! iO.CC 11. oa en c •I"* en «o s- < E C r ^r IT) E » •— • C\J 1 Q C/l LU < < > 1 iJ 00 t\> ^r LU OH ZD CJ3 i—i c> u_ r.V C> 5.03 7.GC 3-0 rRRiSMISSIOH LOSS IN 2E 3.0C a. oa :j c I •a: o lo 3-00 J.J :R9?.SMIS3I3N L?35 1H 2S 11. Qi 162 P 3 •3 « 3 C C7> S- > E c £: Q 8 £ UJ F2 £ I UJ ex: :=> CD 7.\ ii. GO 3- CCi ItfiQ il.QG . 1 163 -J J— _l XJ h- 0) C7> 4-> 0) «=c JC re o E 4-i o c =3 (/> a> a: c a: a: D_ I If) C3 qp *SS01 NOISS!WSN\fbil 164 O _l 1— _J +J 1— c Ol +J s- c 3 O 3 DC DC CO LU DC CD 8P 4SS01 NOISSiWSNVdi 165 to to to - o to _l O' _l c o • T™ to to c U) «d •r- s- E h- to c ■o (O 0) S- CD t— fO -t-J 0> C > (U <: s- 0) E jr c to o +J o LO «o c • o 1— 1 CVJ • /■ o -? O o o o •8 -a c =3 o c C\J «ft e UJ 3 o Q£ z cC < i— i OK h- i o ca < iZ a: UJ o a: o ZD O -* o 8P 'SS01 NOISSIWSNVdl 166 £222mmmh O CM o CO o o E c ouf oo C0 2 < o o "3- to 3 cu n £ 3 a. i to LU a: O QP'SSOI N0ISSIWSNVU1 167 o o o o o <0 CO CO O) o a: o «tf- ~ e c c =3 m a: UJ O 3 3 s » - s £ - r KJ • ^~ BnTTflM 2 - : 1 FIGURE 57 - Antigua to Grand Banks - Bilinear Profil< 169 o o to CJ -J ■a h- - o csJ O o u o z < a: o .o o o 8 S- CD O I0» T o o T o o CM O to gj> *SS0"l N0ISSIWSNVU1 170 o . o CM o o CM Crt 3 to o 0) o cc \Q >> e CD j*. O • r— m Id O CO z < +J E •r" K , — O O O M O zn cr> 00 r o 0> 1^ o o -r o nr o CM o s m LU ce: !3 O »— 1 O ro QP 'SSCH N0ISSIINSNVU1 171 CO CO o CO CO •r- E CO E re s- N CO CD o VO cm ap *SS0"I NOISSIWSNVb'l 172 o o 10 O CJ o o \Q 3 S JC cd m UJ s- e> «o z a: >> CD O o +J N CVJ 10 UJ a: o ID -o O IO i— i I o o I r H^ in o IO ( 2> o — c M 8P 4SS01 NOISSIWSNVdl 174 £ o Z < q: O O O to co o e o CO CO CO c N a: o o o CO CO ZD CD c o o m o in 175 INITIAL DISTRIBUTION LIST No. Copies 1. Defense Documentation Center 2 Cameron Station Alexandria, Virginia 22314 2. Library, Code 0142 2 Naval Postgraduate School Monterey, California 93940 3. Professor Glenn H. Jung, Code 68Jg 1 Department of Oceanography Naval Postgraduate School Monterey, California 93940 4. Assoc. Prof. Alan B. Coppens, Code 61Cz 1 Department of Physics and Chemistry Naval Postgraduate School Mon te rey , Ca 1 i f o rn i a 93940 6. Department of Oceanography, Code 68 3 Naval Postgraduate School Monterey, California 93940 7. Asst. Prof. Robert H. Bourke, Code 68Bf 1 Department of Oceanography Naval Postgraduate School Monterey, California 93940 8. Assoc. Prof. Edward B. Thornton, Code 58Tm 1 Department of Oceanography Naval Postgraduate School Monterey, California 93940 9. Asst. Prof. Stevens P. Tucker, Code 68Tx 1 Department of Oceanography Naval Postgraduate School Monterey, California 93940 10. Acoustic Environmental Support Detachment 1 Office of Naval Research Arlington, Virginia 22217 11. Commanding Officer Fleet Numerical Weather Center Monterey, California 93940 12. Commanding Officer Environmental Prediction Research Facility Naval Postgraduate School Monterey, California 93940 13. Office of Naval Research Code 480 Arlington, Virginia 22217 14. Library, Code 3330 Naval Oceanographi c Office Washington, D. C. 20373 15. Dir. Robert E. Stevenson Scientific Liaison Office, 0NR Scripps Institution of Oceanography La Jolla, California 92037 16. SI0 Library University of California, San Diego P. 0. Box 2367 La Jolla, California 92037 17. Department of Oceanography Library University of Washington Seattle, Washington 98105 18. Commander Naval Weather Service Command Washington Navy Yard Washington, D.C. 20390 19. Oceanographer of the Navy Hoffman II 200 Stoval 1 Street Alexandria, Virginia 22332 20. R. M. Fitzgerald Naval Research Lab, Code 8I76 Washington, D. C. 20375 21. Commander Oceanographic Systems, Atlantic Box 100 Norfolk, Virginia 23511 22. Commander Oceanographic Systems Pacific Box 1390 FPO San Francisco 96610 23. Commanding Officer U. S. Naval Facility FPO New York 09552 24. Commanding Officer U. S. Naval Facility FPO New York 09597 25. LCDR Walter P. Donnelly, USN USS Fox, CG-33 FPO San Francisoc 96601 26. LCDR Kirk E. Evans, USN 1161 Old Springfield Rd. Vandal i a, Ohio 45377 Thesis E766 c.l 1G5S1S Evans A normal mode mode] for estimating low- frequency acoustic^ transmission losV in the deep pee'an. thesE766 A normal mode model for estimating low-f 3 2768 001 89189 8 DUDLEY KNOX LIBRARY