Coastal-trapped waves with stratification, topography, mean flow and bottom friction with complex frequency in Matlab

This set of Matlab mfiles (all having names that begin with “bigc”) can be used to calculate coastaltrapped 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).


The Problem:
We seek to find the free wave solutions for linearized coastal-trapped waves along a straight coast. The depth h(x) is uniform in the alongshore (y) direction. The offshore direction is x (x = 0 is the coast) and the vertical direction is z (z = -h(x) is the bottom), where the origin is at the surface at the coast. Dissipative effects are limited to the infinitesmally thin bottom boundary.
In all cases, subscripted independent variables represent partial differentiation. The mean velocity v0(x, z) is in geostrophic balance with the mean pressure p0 and density ρ0 + ρ1. Turbulent stresses in the x, y directions are given by τ x , τ y , and the bottom stress ( τ0 x , τ0 y ) is related to the interior bottom velocity by τ0 x = ρr uB, τ0 y = ρ r vB, (2a, b) where r(x) 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 (3a) The problem is solved by assuming that where the frequency ω can be complex and the alongshore wavenumber l is real. Using this assumption, the problem can be reduced to a single partial differential equation for pressure (Brink, 2006): This equation and the notation are similar to those of Mooers (1975 a,b), except that these are for hydrostatic conditions.
As noted by Brink, 1982or Dale et al. (2001, this equation admits a spurious solution at ω' = f when v0z = v0x = 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(x) ). In terms of pressure,

x = 0:
There is a choice of two boundary conditions at the coast.

Closed boundary condition:
When there is bottom friction ( UE ≠ 0 ), this condition really only makes sense if the vertical variations in u are negligible over the depth h(0). With no bottom friction, this boundary condition is sensible regardless of the coastal wall height.
Stated in terms of pressure, this condition is: "Open" boundary condition: To mimic an open boundary, we use which does a reasonable job for linear coastal-trapped waves at a distant "open" offshore boundary, where variations should have decayed (see Brink, 1982). To help make sure this is used reasonably, use of this condition requires that the bottom be flat at x = 0.
Stated in terms of pressure, this condition is: This boundary condition can only be expected to work well at subinertial frequencies. At higher frequencies, it can lead to internal wave reflection, rather than the free exit of internal waves.

x = L:
The offshore boundary condition is either the same as the coastal "open" condition (15-16) or the coastal "closed condition" (12)(13)(14). The open condition seems to work well as long as the grid boundary (x = 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 is: or, in terms of pressure: The closed condition is: This form requires that when r ≠ 0, the water is shallow enough that u is essentially depthindependent. I.e., it should typically not be much deeper than about 30 m.
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 to find a complex frequency corresponding to a given wavenumber l. 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, "bigcomp", requires an input array that can be generated by "bigcsetup" and modified by "bigcfinch". When you run "bigcsetup" or "bigcfinch", 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 "bigcomp", 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 symmetric instability is possible, 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 alongshore wavenumber is displayed, and the program displays complex frequency and the inverse resonance parameter (rrr = inverse of the integral of |p| 2 over x 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 is announced at this point. Hit any key to continue, and you will then also get section plots of u, v, w, ρ, energy diagnostics (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, it will then display the next wavenumber and go on from there.
The pressure mode is normalized so that And all of the other (u, v, …) modal structures are consistent with this normalization. This normalization is really only strictly valid for long waves with r = v0 = 0, but it is convenient.
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.

Some Notes on using the software:
The x = 0 and x = L boundary conditions can either be closed (as at a coast), or open. If it is closed, then the water there ought to be shallow if bottom friction is nonzero (with no dissipation, this does not matter). The reason for this is that model balances the bottom Ekman transport in a depth-uniform way at the coastal wall. This makes sense if the water is shallow enough that natural depth variations in u are negligible (hence the model is not imposing any physics that is not there already). This means that with bottom friction, the boundary water depth should be probably less than 50 m, but this depends considerably on conditions, such as the alongshore wavelength (see also Mitchum and Clarke, 1986). If there is no bottom friction, then the coastal wall can have any depth.
If a lateral boundary is open, the bottom should be flat there. In the case of an open boundary condition, the water does not have to be shallow. Care should be taken to make sure that the onshore and offshore boundaries are located far from the main modal structure's location of maximum value.
The numerical grid extends one point outside the x, z domain at each end (nn is the number of x grid points, and mm the number of θ 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 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 hxΔx/(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 alongshore 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 the thermal wind equation. The density field can be left undisturbed at either the offshore or onshore end of the domain.
There is substantial documentation within the mfile "bicomp.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 "bigcuvwr.m" and "bigcomp.m".
All outputs are of arbitrary amplitude (with pressure made real at the top at the coast). Units are consistent within c.g.s., i.e. cm/sec, dyne/cm 2 , gm/cm 3 .
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. The code will virtually always converge to a complex frequency, even if the true result is a real frequency. You need to judge what imaginary value is small enough to be effectively zero, but this is usually pretty obvious. This means, also, that there can be very small energy exchanges in the energy diagnostics when the flow is stable. 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 bigcfinch.m), I find that 0.0001 works out fairly well. This means accuracy to 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).
The input file formats differ slightly between the "bigr*.m" system (for real frequencies only), and this more general "bigc*.m" set. To make the conversions back and forth conveniently, use "c2rinpcon.m" and "r2cinpcon.m".
When the Ertel potential vorticity q = -N 2 f* changes sign relative to f, symmetric instability is possible (Hoskins,1974). So, the routine "bigcrhocal2.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 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 xmax = 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 xmax ≤ 40 km yields reasonable results. This problem arises because the arbitrary forcing used to create the resonance is applied near xmax/2. Trouble can arise when the desired wave has negligible amplitude at this distance offshore.