% auxiliary material describing methods used in GRL paper #2005GL022834
% Thomson et al., "Reflection and tunneling of ocean waves observed at a submarine canyon"

% a typeset pdf of this material is available at 
%	http://science.whoi.edu/users/pvlab/REFLECTION_auxiliary/2005GL022834-methods.pdf
% or by contacting jthomson@whoi.edu and elgar@whoi.edu

% run "latex 2005GL022834-methods.tex" at a unix or linux command prompt 
% to generate a typeset dvi file
% (and then run "dvips 2005GL022834-methods.dvi" to print the file)

\documentclass{article}

% package to include a figure at the end
\usepackage[dvips]{graphicx}

% define reasonable page margins and spacing
\topmargin 0.0cm
\oddsidemargin 0.2cm
\textwidth 16cm
\textheight 21cm
\footskip 1.0cm

% use prefix A on all page, equation, and figure numbers
\renewcommand{\thefigure}{A\arabic{figure}} 
\renewcommand{\theequation}{A\arabic{equation}} 
\renewcommand{\thepage}{A\arabic{page}} 

% footnotes as symbols, not numbers
\renewcommand{\thefootnote}{\fnsymbol{footnote}}


\begin{document}

\noindent{\Large \bf A method to observe the reflection and tunneling\\ 
of ocean waves at a submarine canyon\footnote{This is auxiliary material 
to {\it Geophysical Research Letters}, {\bf XX}, DOI: 10.1029/2005GL022834, 
and is available online at ftp://agu.org/apend/gl/2005GL022834/.  
Hard copies are available by e-mailing jthomson@whoi.edu or elgar@whoi.edu.} \\}

\noindent Jim Thomson\\
{\small WHOI-MIT Joint Program in Oceanography, Woods Hole, Massachusetts, USA\\}

\noindent Steve Elgar\\
{\small Woods Hole Oceanographic Institution, Woods Hole, Massachusetts, USA\\}

\noindent T.H.C. Herbers\\
{\small Naval Postgraduate School, Monterey, California, USA\\}

\baselineskip18pt
\bigskip
 
In contrast to methods used to estimate the amount
of reflection from impermeable structures [{\it Dickson et al.,} 1995],  
estimation of the reflection of a canyon must also account for transmission. 
Here, inverse techniques [{\it Coleman and Li,} 1996] are used to determine the  
reflection ($R_{\rm s}$, $R_{\rm n}$) and transmission ($T_{\rm s}$, $T_{\rm n}$) coefficients,  
as well as the incident directional spectra ($D_{\rm s}$, $D_{\rm n}$)  
and reflected phases ($\psi_{\rm s}$, $\psi_{\rm n}$) 
of the wave fields on the south (s) and north (n) sides of the canyon 
that are most consistent with observations of pressure and velocity. 


The random wave fields on the north (n) and south (s) sides of the canyon consist 
of incident, reflected, and transmitted (from the other side of the canyon) waves 
such that the surface elevations $\eta_{\rm s}$ and $\eta_{\rm n}$ can 
be written as integrals over wave components at each frequency and direction, given by
Eqs. 4 and 5, and repeated here for convenience:  
\begin{equation} 
\eta_{\rm s}  =  \int_{\omega}\!\int_{\theta=-\frac{\pi}{2}}^{\frac{\pi}{2}}  
d_{\rm s} \left( {\rm e}^{i(mx+ly-\omega t)} 
    +\!R_{\rm s} {\rm e}^{i(mx-ly-\omega t + \psi_{\rm s})} \right) 
+  d_{\rm n} \left( T_{\rm n} {\rm e}^{i(mx-ly-\omega t)} \right), 
\end{equation} 
\begin{equation} 
\eta_{\rm n}  =  \int_{\omega}\!\int_{\theta=-\frac{\pi}{2}}^{\frac{\pi}{2}}  
d_{\rm n} \left( {\rm e}^{i(mx-ly-\omega t)}
    +\!R_{\rm n} {\rm e}^{i(mx+ly-\omega t + \psi_{\rm n})} \right)  
+ d_{\rm s} \left( T_{\rm s} {\rm e}^{i(mx+ly-\omega t)} \right) , 
\end{equation} 
\noindent where $d_{\rm s}$ and $d_{\rm n}$ are the complex-valued
differential amplitudes of the incident wave components at 
radian frequency $\omega$ and direction $\theta$
relative to the cross-canyon coordinate $y$. 
The variables $R_{\rm s}$, $R_{\rm n}$ and $T_{\rm s}$, $T_{\rm n}$    
are reflection and transmission coefficients,
$x$ is the along-canyon coordinate,
and $t$ is time.
The along-canyon ($x$) and cross-canyon ($y$)
components of the wavenumber vector are $m=k \sin \theta$ and $l=k \cos \theta$, respectively.
In this two-quadrant ($-\frac{\pi}{2} < \theta < \frac{\pi}{2}$) system, 
the direction of $y$ propagation is given explicitly
by the sign of the exponent in ${\rm e}^{\pm ly}$ because $l$ is always positive.


Assuming specular reflection, the reflected wave fields have amplitudes
$R_{\rm s}d_{\rm s}$ and $R_{\rm n}d_{\rm n}$
that propagate away from the canyon (i.e., ${\rm e}^{\mp ly}$)
with phase shifts $\psi_{\rm s}$ and $\psi_{\rm n}$ relative to the incident wave field.
Conserving energy, the remaining (i.e., nonreflected) portion of the wave field
is transmitted across the canyon with amplitudes
$T_{\rm s}d_{\rm s}$ and $T_{\rm n}d_{\rm n}$.

 
The cross-canyon velocity $v$, along-canyon velocity $u$, and 
pressure $p$ fields induced by surface waves (Eqs. A1 and A2) 
can be determined from a linear, hydrostatic momentum balance
(appropriate for infragravity waves in 20-m water depth)
given by 
\begin{equation}  
\frac{\partial v}{\partial t} = - g \frac{\partial \eta}{\partial y}, 
\quad \frac{\partial u}{\partial t} = - g \frac{\partial \eta}{\partial x}, 
\quad p = \rho g \eta, 
\end{equation} 
where $g$ is gravitational acceleration and $\rho$ is density. 

 
Substituting $\eta_{\rm s}$ and $\eta_{\rm n}$ into the momentum balance (Eq. A3),
using the Fourier transformed result, 
and applying the identities $2 {\bf i} \sin ({\bf i} \alpha) = {\rm e}^{-\alpha} - {\rm e}^{\alpha}$  
and $2 \cos ({\bf i} \alpha) = {\rm e}^{\alpha} + {\rm e}^{-\alpha}$
yields the following expressions for the frequency cross-spectra 
of the colocated pressure and velocity time series south (s) of the canyon  
\begin{equation}  
\langle p_{\rm s}(\omega)^{ } \!\cdot\! u_{\rm s}^{\star}(\omega) \rangle\!=\!\rho g \kappa\!\int_{\theta} \sin \theta \left[ D_{\rm n} T_{\rm n}^2 
+ D_{\rm s} \left( 1+R_{\rm s}^2+2R_{\rm s}\cos \psi_{\rm s} \right) \right] d \theta, 
\end{equation} 
\begin{equation} 
\langle p_{\rm s}^{ }(\omega) \!\cdot\! v_{\rm s}^{\star}(\omega) \rangle\!=\!\rho g \kappa\!\int_{\theta} \cos \theta \left[ -D_{\rm n} T_{\rm n}^2 
+ D_{\rm s} \left(1 - R_{\rm s}^2  + {\bf i}2R_{\rm s} \sin \psi_{\rm s} \right) \right] d \theta, 
\end{equation} 
\begin{equation} 
\langle u_{\rm s}^{ }(\omega)\!\cdot\! v_{\rm s}^{\star}(\omega) \rangle\!=\!\kappa^2\!\int_{\theta} \sin\theta\cos\theta \left[ -D_{\rm n} T_{\rm n}^2 
+ D_{\rm s} \left(1\!-\!R_{\rm s}^2\!+\!{\bf i} 2R_{\rm s} \sin \psi_{\rm s} \right) \right] d \theta , 
\end{equation}
where $\kappa = \frac{gk}{\omega}$, $\star$ is the complex conjugate, 
and  $\langle \rangle$ is the expected value.
Similarly, the auto-spectra are
\begin{equation}  
\langle p_{\rm s}^{ }(\omega) \!\cdot\! p_{\rm s}^{\star}(\omega) \rangle\!=\!(\rho g)^2\!\int_{\theta} \left[ D_{\rm n} T_{\rm n}^2 
+ D_{\rm s} \left( 1+R_{\rm s}^2+2R_{\rm s}\cos \psi_{\rm s} \right) \right] d \theta, 
\end{equation} 
\begin{equation} 
\langle u_{\rm s}^{ }(\omega) \!\cdot\! u_{\rm s}^{\star}(\omega) \rangle\!=\!\kappa^2\!\int_{\theta} \sin^2\theta \left[ D_{\rm n} T_{\rm n}^2 
+ D_{\rm s} \left(1 + R_{\rm s}^2  + 2R_{\rm s} \cos \psi_{\rm s} \right) \right] d \theta, 
\end{equation} 
\begin{equation} 
\langle v_{\rm s}^{ }(\omega) \!\cdot\! v_{\rm s}^{\star}(\omega) \rangle\!=\!\kappa^2\!\int_{\theta} \cos^2\theta \left[ D_{\rm n} T_{\rm n}^2 
+ D_{\rm s} \left(1 + R_{\rm s}^2 - 2R_{\rm s} \sin \psi_{\rm s} \right) \right] d \theta ,  
\end{equation}
The incident wave directional spectra $D_s$ and $D_n$ are defined as  
\begin{equation} 
D_s(\omega,\theta) = \frac{\langle d_s \cdot d_s^{\star} \rangle }{d \omega d \theta},
\end{equation}
\begin{equation}
D_n(\omega,\theta) = \frac{\langle d_n \cdot d_n^{\star} \rangle }{d \omega d \theta}. 
\end{equation} 
The incident waves fields on the north and south sides
of the canyon are independent of each other, and thus 
$\langle d_s \cdot d_n^{\star} \rangle=
\langle d_n \cdot d_s^{\star} \rangle=0.$

Expressions for cross- and auto-spectra at the north side of the canyon are 
obtained by exchanging all subscripts (s$\leftrightarrow$n) in the expressions  
above (Eqs. A4-A9) and changing the sign of the integrand in Eqs. A5 and A6. 

 
The real-valued terms describe the progressive wave field, and the 
imaginary (${\bf i}$) terms  
(i.e., the quadrature in the cross-spectra $\langle p \cdot v^{\star} \rangle$ 
and $\langle u \cdot v^{\star} \rangle $)  
describe the partial standing wave patterns 
owing to sums of incident and reflected waves.  
In practice, the expressions for cross- and auto-spectra apply only 
to observations near the canyon walls,  
because over large distances (i.e., many wavelengths)  
standing wave patterns are obscured within a finite-width frequency band. 

  
For computational efficiency, the incident directional spectra (Eqs. A10 and A11) at 
each frequency band are modeled as [{\it Donelan et al.,} 1985]
\begin{equation} 
D_{\rm s}(\theta) = M_{\rm s} \cos^{2S_{\rm s}} \left( \frac{\theta-\Theta_{\rm s}}{2} \right),
\end{equation}
\begin{equation}
D_{\rm n}(\theta) = M_{\rm n} \cos^{2S_{\rm n}} \left( \frac{\theta-\Theta_{\rm n}}{2} \right),
\end{equation} 
where $\Theta_{\rm s}$ and $\Theta_{\rm n}$ are the centroidal directions,
$S_{\rm s}$ and $S_{\rm n}$ describe the spread in direction about the centroid,
and $M_{\rm s}$ and $M_{\rm n}$ are the spectral peak values.
The results are insensitive to the specific unimodal shape used for the incident directional spectra. 
The centroidal directions $\Theta_{\rm s}$ and $\Theta_{\rm n}$  
were used to separate the data sets (at each frequency band) into normally ($|\Theta| < 20^\circ$)
and obliquely ($|\Theta| > 30^\circ$) incident waves.
Assuming directionally narrow spectra, 
the reflection ($R_{\rm s}, R_{\rm n}$) and 
transmission ($T_{\rm s}, T_{\rm n}$) coefficients are
assumed to be independent of direction at each frequency.
The spreading parameters ($S_{\rm s}$, $S_{\rm n}$)  
were used to weight individual estimates of $R^2$ before averaging 
over the collection of data runs (e.g., the symbols in Figure 2).


The phase shifts $\psi_{\rm s}$ and $\psi_{\rm n}$
of the reflected waves relative to the incident waves were
allowed to vary over each directional spectrum by assuming
\begin{equation}
\psi_{\rm s} (\theta) =2 \Delta y_{\rm s} k \cos \theta,
\end{equation}
\begin{equation}
\psi_{\rm n} (\theta) =2 \Delta y_{\rm n} k \cos \theta,
\end{equation}
where $\Delta y_{\rm s}$ and $\Delta y_{\rm n}$ are the (unknown) distances
between the reflector and the instrument locations (i.e., $\psi$ is the
phase change associated with propagating toward and back from the reflector).
The resulting inverse estimates of $\Delta y_{\rm s}$ and $\Delta y_{\rm n}$
are consistent with the measured distances (Figure A1).  
The small deviations from theory may be caused by incorrectly assuming all the reflection 
takes place at a single location.  

 
Assuming energy is conserved (i.e., $R_{\rm s}^2+T_{\rm s}^2=1$, $R_{\rm n}^2+T_{\rm n}^2=1$),  
the inverse method finds the north (n) and south (s) values of  
$\Theta, S, M, R,$ and $\Delta y$ that are most consistent with the  
cross- and auto-spectra of the observed time series 
by minimizing a normalized root-mean-square error [{\it Dickson et al.,} 1995]   
\begin{equation}                                            
\epsilon = \sqrt{\frac{\sum(obs-derived) \cdot (obs-derived)^{\star}} 
{\sum(obs)(obs)^{\star}}},
\end{equation}                                    
where $\sum$ is the sum over the six spectral values from the south side (Eqs. A4-A9) and 
the six spectral values from the north side.
Applying recent improvements to Newton's method [{\it Coleman and Li,} 1996], 
the inverse algorithm begins with an initial guess for each unknown, 
and solves a locally linearized version of the equations for the cross- and auto-spectra
to find the small change in each unknown that produces the greatest reduction in $\epsilon$
(i.e., the method iterates down-slope in $\epsilon$ until the minimum is found).   
Initial guesses for $\Theta$, $S$, and $M$ are provided by estimates of directional moments of the
wave field based on the measured pressure and velocities [{\it Kuik et al.,} 1988],
and the initial guesses for $R$ are based on long-wave theory [{\it Kirby and Dalrymple,} 1983].
The results are not sensitive to the initial values, and the same inverse solutions 
are obtained with random initial guesses
(although computational time is increased). 


The residual errors (Eq. A16) approximately follow the expected $\chi^2$ distribution  
of fluctuations in spectral estimates from finite length data records [{\it Wunsch,} 1996]. 
The percent of observed variance captured by the inverse method (defined as $100 \times [1-\epsilon^2]$) 
is 90\% when averaged over all infragravity frequency bands for all 50 data sets. 



\newpage
\baselineskip12pt
\section*{References}
\begin{description}
\item[] Coleman, T., and Y. Li (1996), An interior, trust region approach for nonlinear minimization subject to bounds, {\it SIAM J. Optimization}, {\it 6}, 418-445. 
\item[] Dickson, W., T.H.C. Herbers, and E. Thornton (1995), Wave reflection from breakwater, {\it J. of Water., Port, Coast., Ocean Eng.} {\bf 121}, 262-268. 
\item[] Donelan, M., J. Hamilton, and W. Hui (1985), Directional spectra of wind-generated waves, {\it Phil. Trans. Roy. Soc. Lon.}  {\bf A315}, 509-562. 
\item[] Kuik, A., G. Van Vledder, and L. Holthuijsen (1988), A method for the routine analysis of pitch-and-roll buoy wave data, {\it J. Phys. Oceanog.} {\bf 18}, 1020-1034. 
\item[] Wunsch, C. (1996), {\it The Ocean Circulation Inverse Problem} (Cambridge University Press, Cambridge. 
\end{description}
 
 
\begin{figure} 
\begin{center}
\noindent \includegraphics{2005GL022834-methods-figureA1} 
\caption{\sf Distance (m) between the reflector and the instrument location versus frequency (mHz)   
and period (s) at the (a) south and (b) north sides of the canyon.  
Symbols are estimates from the nonlinear inverse method and  
solid lines are the measured distances from the steepest portion of  
each canyon slope to the instrument site at that side.
Instrument locations were determined within $\pm$ 10 m
(the width of the solid lines) with differential GPS. 
The theory assumes that waves propagate from the  
instrument site to a vertical-canyon-wall reflector along a line of constant $y$ and back,  
neglecting possible phase shifts at the sloped walls.   
Vertical lines are $\pm$ one standard deviation of the estimates.}  
\end{center}
\end{figure} 

\end{document}
