Seamount-or Lake / Basin-trapped waves with stratification , topography , mean flow and bottom friction in Matlab

This set of Matlab mfiles (all having names that begin with “bigs”) can be used to calculate seamount-trapped (or basin-trapped) wave modal structures and dispersion curves under very general circumstances. A complex frequency is allowed, so that instability and damping can be accounted for directly. Modal structures and energy diagnostics are provided. For most applications, the code is only useful for subinertial wave frequencies (i.e., the real part of wave frequency is smaller than the Coriolis parameter). For interpreting the model results, see Brink (1989), which deals with the case with no mean flow or finite bottom friction. The present code was developed independently of the Fortran code used in that publication.


The Problem:
We seek to find the free wave solutions for linearized seamount-trapped waves at a circular seamount.The depth h(r) is uniform in the azimuthal (θ) direction.The outward (radial) direction is r and the vertical direction is z.The bottom is at z = -h(r) and the origin is at the seamount/basin center.Dissipative effects are limited to the infinitesmally thin bottom boundary layer.In a seamount configuration, the shallowest water is at the center and the outer boundary is open (and the bottom is flat there).In a lake configuration, the deepest water is at the center, and the outer boundary is closed.In a basin case, the deepest water is at the center, but the outer boundary is open and the bottom is flat there.
The mean azimuthal flow v0(r, z) obeys -v0 (f + v0 / r) = -ρ0 -1 p0r (2a) so that v0z (f + 2v0 / r) = -g ρ0 -1 ρ1r . (2b) Turbulent stresses in the r, θ directions are given by τ r , τ θ , and the bottom stress ( τB r , τB θ ) is related to the interior bottom velocity by τB r = ρ0 R uB, τB θ = ρ0 R vB, (3 a, b) where R(r) is a bottom resistance coefficient having units of velocity, and a subscript "B" denotes a quantity evaluated at the bottom.We also define N 2 = -gρ0 -1 ρ1z and (4a) The problem is solved by assuming that p = p´ (r, θ) exp[i(ωt + mθ)], etc., where the frequency ω can be complex and the azimuthal wavenumber m is a real integer.Henceforth, the "primes" will be dropped.Using this form, the problem can be reduced to a single partial differential equation for pressure (analogous to Brink, 2006) 0 = ω´ prr + prz (-2 sω´2) + pzz ω´ (s 2 + G/N 2 ) (6) where This equation and the notation are similar to those of Mooers (1975 a, b), except that these are in cylindrical coordinates and apply to hydrostatic conditions.
The mean azimuthal velocity field is taken to be a modified Gaussian in r and z, centered at r= rv, z = -zv.
where <v> is the peak value, and the horizontal scale D r has different values onshore and offshore of rv, while the vertical scale D z is allowed to be different upward and downward of zv.
The function F(r) assures that v0 → 0 as r → 0. Thus, if you chose a large D r for r < rv, the mean flow will simply be in solid-body rotation for r < rv.If you desire to have a different initial velocity form, you can simply change the code in file bigsvzcal2.m,around lines 35-81.Velocity gradients and density are computed numerically given v0, so that no other code need be changed.
As noted by Brink, (1982) or Dale et al. (2001), equation ( 6) admits a spurious solution at ω' = f when v0 = 0.The presence of this mode also distorts modal structures and dispersion curves for nearby frequencies.Typically, this spurious pressure mode is vertically uniform and it decays only slowly with distance offshore.
The seamount problem involves having minimum water depth at the center, and then increasing depth to an outer, flat region.An open outer boundary condition would apply.The real part of frequency would be > 0 in the northern hemisphere in the absence of a mean flow.The basin/lake problem would have maximum depth at the center and shoaling outward.Either a closed (lake) or open (basin) boundary condition would apply at the outer boundary.The real part of frequency would be < 0 in the absence of a mean flow.
The boundary conditions are as follows: where δ = 1 for a free surface condition and δ = 0 for a rigid lid condition.In terms of pressure, where the Ekman transport is given as and the subscript "B" denotes a quantity evaluated at the bottom (along z = -h(r) ).In terms of pressure, The boundary condition at the center (r = 0) is simply

Open outer boundary condition (iobc = 1)
The outer boundary condition is an adaptation of the approximate open condition introduced by Brink (1982): This condition seems to work well as long as the outer grid boundary (r = L) occurs far offshore of where the modal solutions have substantial structure.The condition, as used, requires that the bottom is flat at the offshore boundary.The condition in terms of pressure is: 0 = -ω´ prr + sω´ pzr + pr (ω´Gr / G -ω´/ r -ω´r -mf1/r) (16)

Closed outer boundary condition:
If there is no bottom friction, this condition is simply This boundary condition is then valid regardless of the water depth at r = L.If there is bottom friction (R ≠ 0), then the condition is where the approximation is only valid if the water is shallow enough that u does not vary substantially over depth h.Thus, it is important that if you have a closed outer boundary condition (as in a lake), the coastal wall is not very high (certainly less than about 30 m).
Expressed in terms of pressure, the closed outer boundary condition in either case becomes The problem is now expressed in stretched coordinates that map the problem into a rectangular domain (e.g.Wang and Mooers 1976).Specifically, The problem is then solved numerically via resonance iteration.An arbitrary forcing is applied through the surface boundary condition (essentially, a localized wind stress curl near the grid center), and ω is searched for a maximum response so as to establish a complex frequency corresponding to a given wavenumber m.This search makes use of the Matlab function "fminsearch.m".

Using the software:
Examples of what your Matlab screen will look like as you run the different mfiles are shown at the end of this document.The main routine, "bigsomp", requires an input array that can be generated by "bigssetup" and modified by "bigsfinch".When you run "bigssetup" or "bigsfinch", you simply need to answer the questions posed.Do not forget to use square brackets when entering more than one number on a line.
When running "bigsomp", you will first be shown a few simple parameter values.You also get a section plot showing topography, stratification (minus a constant background value), and the mean flow.If f*/f is negative anywhere, f* will also be plotted (in red) on the density plot.If you are allowing pauses during the run, you are advised of pauses, and need to hit a key each time to continue.During iteration toward a solution, the azimuthal wavenumber m is displayed, and the program displays complex frequency and the inverse resonance parameter (rrr = inverse of the integral of r|p| 2 over r and z) for each ω estimate.When you are at a resonance (free mode), rrr becomes very small.You can thus watch the convergence onto a complex free modal frequency.
Convergence to your given tolerance is then announced, and the free mode pressure is the plotted as a section.A pause can also be announced at this point.Hit any key to continue, and you will then also get section plots of u, v, w, ρ, energy diagnostics in normalized units (energy pool sizes and transfers: the dissipation integral should be a negative quantity), and line plots of p, v, u at the surface.If you set the program to run more than one frequency/wavenumber, it will then display the next wavenumber and go on from there.
Before any plots are made, the pressure field is normalized (following Brink, 1999) so that where η = 1 for a closed outer boundary (lake) and η = 0 for an open outer boundary (seamount) (This is controlled by your choice of iobc).This means that, in the following plots, pressure has units of (cm -½ ), velocity has units of (cm -3/2 sec gm -1 ), and density has units of (cm -5/2 sec 2 ).The energy diagnostics will also have consistent peculiar units.This form is not technically valid when either R ≠ 0 and/or <v> ≠ 0, but it is applied nonetheless for consistency.
If you have set the program to compute more than one frequency/wavenumber pair, the program will also draw a complex dispersion curve once all pairs have been obtained.In this case, you will be asked at the end of the run whether you want to save the dispersion curve as a .matfile.

Some Notes on using the software:
At the outer boundary, the bottom should be flat if you are using an open boundary condition, and the water can be any depth there.In the seamount case, care should be taken to make sure that the outer boundary is located far from the topography and the modal structure's location of maximum amplitude.In the case of a closed outer boundary (lake), the bottom need not be flat at the boundary.The numerical grid extends one point outside the r, z domain at each end (nn is the number of r grid points, and mm the number of z (or, equivalently φ) grid points.There is also a grid point beyond each boundary.Thus, for example, nn = 5 gives one interior grid point, 2 on boundaries and 2 outside the domain.
In choosing a grid size, caution is advised.One should experiment with grid sizes in order to be confident of a good answer.I find that nn = 80-120, mm = 20-30 seems to work well in many coastal cases.Because seamount-trapped waves tend to be bottom-trapped, you may want to use better vertical resolution than this, like mm = 40-60.The number of operations or amount of memory needed goes something like nn*mm*mm.
There is a test for hydrostatic consistency (the ratio computed is |hrΔr/(hΔφ)| ), and results seem best if the ratio given is less than around 1-5. .Further, steep bottom topography can cause problems.It appears that a reasonable criterion is to pick nn such that The mean azimuthal flow is assumed to be Gaussian with a linear correction inshore of the core (with different onshore, offshore, up and down length scales, as well as a variable center position and depth).This provides a fair bit of flexibility.Once the velocity field is computed, the density structure is computed using a variant of the thermal wind equation ( 2b).The density field can be left undisturbed at either the outer or inner end of the domain.The program will alert you (by blue "+" signs) to any regions where your choices have led to unstable vertical density gradients, i.e., when (ρ0 + ρ1)z > 0.
There is substantial documentation within the mfile "bigsomp.m"for further detail.
Results for pressure, velocity and density are plotted.The smaller of the real or imaginary is suppressed if it is too small (maximum absolute value less than about 1000 times less than the other's maximum absolute value, as it is set now).This ratio is controlled by the value of "ratnplot" in "bigsuvwr.m"and "bigsomp.m".
The search can be very inefficient if you give it 0 as the starting guess for real or imaginary frequency.Better to give it a very small nonzero value.For example, if you are seeking a mode with a real frequency (say, when R and v0 are both zero), then a good starting guess might be 1×10 -5 + i 1×10 -13 sec -1 rather than simply 1×10 -5 + i 0 sec -1 .Similarly, if you are seeking a damped wave frequency (e.g., if R ≠ 0 and v0 = 0), then it would be well to start with a guess like (1 + i)×2×10 -5 sec -1 , rather than (2×10 -5 + i 1×10 -10 ) sec -1 , for example (i.e., use similar real and imaginary magnitudes, rather than a large and very small value).Similarly, if you are seeking an unstable mode, you might want to start searching in a similar way, e.g., (1 -i)×5×10 -6 sec -1 .The magnitude of these numbers are just for illustration.The relative values of the real and imaginary parts are the real point.The search will usually converge on some value, but it may not be meaningful.The search might also miss meaningful solutions if the initial guess is too far off.
In considering the outputs, you need to judge whether the imaginary value is small enough to be effectively zero, but this is usually pretty obvious.Looking at the energy diagnostics can often tell you if an answer is silly.The nonzero imaginary frequency means, also, that there can be very small energy exchanges in the energy diagnostics when the flow is stable.Further, if the contoured solutions show a concentration of contours near-surface in the center of the grid (where the arbitrary forcing occurs), then the result is likely bogus.Again, you need to decide what is effectively zero.Also, for the seamount case, the modal structure for p can sometimes look quite different from that for v, for example.This should not cause any problems.
By the same token, if you give the algorithm an initial estimate for real or complex frequency that has few places of accuracy (say 1e-5), the Matlab search code will look around at fairly large increments (perhaps ± 20%), and the search can miss nearby solutions.If you provide more places of accuracy (say 1.01e-5), the search will be initially confined to more nearby locations (perhaps ± 3%), and it becomes much more likely that nearby solutions will be found.
The search can sometime miss a complex frequency, for example if searching all takes place only near the real axis, but the true result has a substantial imaginary part.If it is a good solution, the inverse resonance parameter in the output (rrr) should be several orders of magnitude lower than neighboring values.Use a smaller accuracy tolerance (e) or start the search again farther from the real axis.
For an accuracy estimate ("e" option in bigsfinch.m),I find that 0.0001 works out fairly well.This means a nominal accuracy of 0.01% of the absolute value of the initial frequency guess.
The bottom stress is assumed to be proportional to the bottom velocity.For the errors associated with this assumption in the presence of a mean flow, see Brink (1997).
If more than one frequency is computed, there is an option to save the "dispersion curve" information (wavenumbers and complex frequencies).
When the Ertel potential vorticity q = -N 2 f* changes sign relative to f, symmetric instability is possible (Hoskins, 1974).So, the routine "bigsrhocal2.m"tests to see if f*/f is negative anywhere.If this happens, f* is contoured in red on the same plot as mean density, and a message appears on your Matlab screen.If f*/f does not become negative, it is not contoured, and no message about it is given.Experience with the program suggests that modes calculated when f* changes sign are of questionable merit.Structures are then often complex and have short scales near the region of negative f*.It seems probable that if, in nature, a flow is symmetrically unstable, it will quickly mix to a neutral state with zero or positive f*/f.So, even a "good" modal solution for a symmetrically unstable case is of questionable value.
Given the presence of the near-inertial spurious mode, and given that the "open" boundary condition is only reasonable for subinertial frequencies, the seamount model results should only be trusted for subinertial frequencies, |Re(ω)| <| f |.If, in computing a dispersion curve, you start to obtain real frequencies greater than the inertial, the superinertial results should be discarded.
Once you have an output showing the modal structure, you may want to decide which radial mode number (n) it is.The easiest way is to look at the plot of pressure along the bottom (top plot in Matlab figure 6).The mode number for a seamount-trapped wave is usually the number of zero crossings in bottom pressure, including the zero at r = 0.For a lake-or basin wave, it is better to count zero crossings in v.This tends to work well when R = v0 = 0.With friction or a mean flow, making this identification can become trickier, but you still should look at the real part of pressure only.
All frequencies are in radians per second, and all wavenumbers in radians per cm.Radians do not have units, so these units could be written as "1/sec" for example.
You can get spurious results if the domain is too large relative to the natural scale of the wave under consideration.For example, with a domain size L = 70 km, trying to find a flat-bottom internal Kelvin wave with internal Rossby radius = 3 km yields bad results.Doing the same calculation with L ≤ 40 km yields reasonable results.This problem arises because the arbitrary forcing used to create the resonance is applied near L/2.Trouble can arise when the desired wave has negligible amplitude at this distance offshore.

Mfiles in this package:
This collection of mfiles includes the following:

bigsomp.m
The driver for the seamount-and basin-trapped wave calculations.Carries out all input functions, calls the functions for computing mean density and velocity fields, and for searching for resonant frequencies.Also does some output and graphics functions.

bigssetup.m
This mfile sets up an input array for "bigsomp", based on user responses to questions.An example of its use is given below.

bigsfinch.m
This mfile changes an existing "bigsomp" input array, without having to generate an all new array.An example of its use is given below.

bigsdep.m
Computes a gridded depth profile, along with its derivatives.

bigsconchk.m
Grid consistency check to see if the bottom slope is large relative to vertical spacing.

bigsvzcal2.m
Calculates the mean velocity field and its vertical derivative, from inputs.

bigsrhocal2.m
Calculates the modifications to the basic density field caused by the mean flow.

bigsdrver.m
Mainly called from fminsearch.This sets up matrix equations that are equivalent to the field differential equation and its boundary conditions.Also solves the equation, and computes an integral of the response to evaluate convergence.

bigssigtoz.m
Converts a field from sigma coordinates to z coordinates.

bigsengdiag.m
Does the energy diagnostics for each solution.

bigsztosig.m
Converts z coordinates to sigma coordinates.

bigscc.m
Computes coefficients for the center boundary condition.

bigsoff.m
Computes coefficients for the outer boundary condition.

bigsbot.m
Computes coefficients for the bottom boundary condition.

bigssbc.m
Computes coefficients for the surface boundary condition.

bigsmatco.m
Computes coefficients for interior points on the finite difference grid.

bigsint.m
Computes the integral of the response magnitude.
Also, a sample input file for Fieberling Guyot (fieb0.mat) is also included.