Island-trapped waves with stratification , topography , mean flow and bottom friction in Matlab

This set of Matlab mfiles (all having names that begin with “bigi”) can be used to calculate island-trapped wave modal structures and dispersion curves under very general circumstances for a circular island. 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 (1999), 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 island-trapped waves at a circular island.The depth h(r) is uniform in the azimuthal (θ) direction.The offshore (radial) direction is r (r = r0 is the coast) and the vertical direction is z.The bottom is at z = -h(r) and the origin is at the island center.Dissipative effects are limited to the infinitesmally thin bottom boundary layer.

(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 = ρ0R 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., (5) 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 that in Brink, 2006): 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 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.
If you desire to have a different initial velocity form, you can simply change the code in file bigivzcal2.m,around lines 35-78.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 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, A closed boundary condition gives no net flux across the coastal boundary: When there is bottom friction (UE ≠ 0 ), this condition really only makes sense if the vertical variations in u are negligible over the depth h0 = h(r0).With no bottom friction, this boundary condition is sensible regardless of the coastal wall height.
Stated in terms of pressure, this condition is: The offshore 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) (17) 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, "bigiomp", requires an input array that can be generated by "bigisetup" and modified by "bigifinch".When you run "bigisetup" or "bigifinch", 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 "bigiomp", 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 arbitrary 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 (Brink, 1999) so that 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 ).Relations in Brink (1999) can be used to estimate wind coupling coefficients when R = 0 and <v> = 0, given his normalized modal form.
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:
The r =r0 boundary condition is closed.This causes no issues when R = 0, since it just imposes u = 0 at each depth, a reasonable condition at a coastal wall.When R ≠ 0 , bottom friction is applied, and the bottom Ekman transport that reaches the coastal wall is balanced by an evenly distributed interior flow.Thus, there is an assumption that the water is shallow enough that |uz /u| is small.This would be expected if (f/N) L r is large relative to h0, where L r is a representative modal scale in the r direction.This means that, with bottom friction, the boundary water depth should be probably substantially less than 50 m.Further, see Mitchum and Clarke (1986) for aspects of boundary layer theory that affect one's choice of wall height.If there is no bottom friction, then the coastal wall can have any depth.
At the offshore boundary, the bottom should be flat, but the water does not need to be shallow Care should be taken to make sure that the offshore boundary is located far from the main modal structure's location of maximum amplitude.
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 and island cases.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 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 offshore or onshore end of the domain.The program will alert you (by blue "+" signs) to any regions where your choices have led to gravitationally unstable vertical density gradients, i.e., when (ρ0 + ρ1)z > 0.
There is substantial documentation within the mfile "bigiomp.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 "bigiuvwr.m"and "bigiomp.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 the numbers here 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.
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 bigifinch.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 "bigirhocal2.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 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 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.
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.
A file is included with the input information for a realistic Bermuda case (berm0.mat).

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

bigiomp.m
The driver for the island-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.An example is given below.

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

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

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

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

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

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

bigidrver.m
Mainly called from fminsearch.This sets up the 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.

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

bigiengdiag.m
Does the energy diagnostics for each solution.

bigiztosig.m
Converts z coordinates to sigma coordinates.

bigicc.m
Computes coefficients for the coastal boundary condition.

bigioff.m
Computes coefficients for the offshore boundary condition.

bigibot.m
Computes coefficients for the bottom boundary condition.

bigisbc.m
Computes coefficients for the surface boundary condition.

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