2.1. Non-dimensional compressible-fluid equations

Atmospheric observations following the Tonga event demonstrate that the pressure disturbances in the atmosphere can travel uninterrupted around the globe for several days (Matoza Reference Matoza2022), suggesting that dissipative effects may be neglected. From this observation the fluid is assumed to obey the non-dimensional compressible Euler equations (i.e. friction and heat losses are neglected),

(2.1)\begin{equation} St\frac{\partial}{\partial t} \left[\begin{matrix} \rho \\ \rho \boldsymbol{v} \\ \rho e_T \end{matrix}\right] + \boldsymbol{\nabla} \boldsymbol{\cdot} \left[\begin{matrix} \rho\boldsymbol{v} \\ \rho\boldsymbol{v}\otimes\boldsymbol{v} + Eu p {\boldsymbol{\mathsf{I}}}_3 \\ \rho\boldsymbol{v} e_T + Eu p\boldsymbol{v} \end{matrix}\right] = \left[\begin{matrix} 0 \\ \rho\boldsymbol{g}/Fr^2 \\ \rho\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{g}/Fr^2 \end{matrix}\right],\end{equation}

where $\rho$, $\boldsymbol {v}$, $p$, $e_T$ are, respectively, the density, velocity vector, thermodynamic pressure and specific total energy of the fluid, with $e_T \equiv \boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {v}/2 + e$, where $e$ is the specific internal energy of the fluid. Here $\boldsymbol {\nabla }\boldsymbol {\cdot }()$ is the divergence operator, $t$ is time and $\boldsymbol {v}\otimes \boldsymbol {v}$ is the velocity dyadic product. Non-dimensional numbers $St\equiv \ell _o /(t_o u_o)$, $Eu\equiv p_o/(\rho _o u_o^2)$, $Fr\equiv u_o/\sqrt {g_o\ell _o}$ are, respectively, the Strouhal, Euler and Froude numbers, where $t_o$, $u_o$, $\ell _o$, $p_o$, $\rho _o$, are reference time, speed, length, pressure, density, scales, respectively, and $g_o$ is a reference gravitational acceleration. Term $\boldsymbol {g}$ is the constant magnitude and direction gravitational acceleration vector non-dimensionalised by $g_o$ (tidal effects are neglected). Term ${\boldsymbol{\mathsf{I}}}_n$ is the identity tensor in $n$ dimensions. Without loss of generality, the reference time, pressure and speed are set to $t_o = \ell _o/u_o$, $p_o = \rho _o u_o^2$ and $u_o = \sqrt {g_o\ell _o}$, resulting in $St = Eu = Fr = 1$.

2.2. Depth averaging

Consider the geocentric spherical coordinate system with origin $O$ at the centre of mass of the Earth, unitary basis $(O;(\boldsymbol {e}_{\mathcal {r}},\boldsymbol {e}_\theta,\boldsymbol {e}_\varphi ))$, coordinates $(\mathcal {r},\theta,\varphi )$ and position vector $\boldsymbol {r}=\mathcal {r}\boldsymbol {e}_\mathcal {r}$, where the radius $\mathcal {r}$ is defined as the distance from the origin, the colatitude $\theta$ is defined as the polar angle measured from the North Pole and the longitude $\varphi$ is defined as the angle in the equatorial plane measured from the prime meridian positive to the east. Integrating (2.1) in the radial direction $\boldsymbol {e}_\mathcal {r}$ between the two surfaces bounding a fluid layer results in the depth-averaged equations. Let $\mathscr {Z}_T(t,\mathcal {r},\theta,\varphi )=0$ and $\mathscr {Z}_B(t,\mathcal {r},\theta,\varphi )=0$ define, respectively, the top and bottom surfaces. Assuming these surfaces to be smooth (i.e. neglecting breaking waves and imposing $\partial \mathscr {Z}_i/\partial \mathcal {r}\neq 0$), the implicit function theorem states that these can be redefined as

(2.2)\begin{equation} \mathscr{Z}_i(t,\mathcal{r},\theta,\varphi)\equiv\eta_i(t,\theta,\varphi)-\mathcal{r}=0, \quad \text{with} \ i\in\{B,T\}. \end{equation}

The evolution of these surfaces represents the kinematic boundary condition given by

(2.3)\begin{equation} \frac{{\rm D}\mathscr{Z}_i}{{\rm D}t} =\frac{\partial\eta_i}{\partial t}+\boldsymbol{v}|_{\eta_i} \boldsymbol{\cdot}\nabla^{\eta_i}\mathscr{Z}_i=0, \end{equation}

where $\boldsymbol {v}=[v_\mathcal {r},v_\theta,v_\varphi ]^{{\rm T}}$ is the velocity vector, the in-plane gradient operator for an arbitrary radius $a$ is defined as

(2.4)\begin{equation} \nabla^a\phi\equiv\left[\frac{\partial\phi}{\partial\mathcal{r}}, \frac{1}{a}\frac{\partial \phi}{\partial \theta}, \frac{1}{a\sin\theta}\frac{\partial \phi}{\partial \varphi}\right]^{{\rm T}}. \end{equation}

Throughout the text $\phi$ and $\boldsymbol {f}$ denote an arbitrary scalar and vector field, respectively.

Let $\langle \phi \rangle$ and $\bar {\phi }$ denote, respectively, the linear and logarithmic depth averages,

(2.5a,b)\begin{equation} \langle \phi \rangle\equiv\frac{\displaystyle \int_{\eta_B}^{\eta_T}\phi \,\mathrm{d} \mathcal{r}}{\mathcal{H}} \quad\text{and}\quad \bar{\phi}\equiv\frac{\displaystyle \int_{z_B}^{z_T}\phi\,\mathrm{d} z}{\mathcal{L}},\end{equation}

where $\mathcal {H}\equiv \eta _T-\eta _B$ is the non-dimensional layer thickness, $\mathcal {L}\equiv \ln (\eta _T/\eta _B)$ and $z=\ln (\mathcal {r})$, $z_B = \ln (\eta _B)$, $z_T = \ln (\eta _T)$.

Appendix B derives a general expression for the depth-averaged divergence of a vector field in spherical coordinates and its approximation under the thin-layer assumption. Using Leibniz integration rule and the results of Appendix B, the depth average of the time derivative and divergence are written as

(2.6)\begin{gather} \int_{\eta_B}^{\eta_T}\frac{\partial \phi}{\partial t}\, \mathrm{d} \mathcal{r} =\frac{\partial }{\partial t}(\mathcal{H}\langle \phi \rangle)- \left[\!\!\left[ {\phi\frac{\partial\eta}{\partial t}} \right]\!\!\right], \end{gather}

(2.7)\begin{gather}\int_{\eta_B}^{\eta_T}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{f}\,\mathrm{d} \mathcal{r}= {\nabla}_{{\perp}}^\mathcal{R} \boldsymbol{\cdot}(\mathcal{L}\mathcal{R}{\bar{\boldsymbol{f}}}_{{\perp}}) -[\![\,{\boldsymbol{f}\boldsymbol{\cdot}\nabla^\eta\mathscr{Z}} ]\!]+ 2\mathcal{L}\bar{\boldsymbol{f}}\boldsymbol{\cdot}\boldsymbol{e}_\mathcal{r}, \end{gather}

where $\mathcal {R}$ is an arbitrary constant radius, the subscript $\perp$ denotes in-plane components of the vector (orthogonal to $\boldsymbol {e}_\mathcal {r}$), defined as ${\boldsymbol {f}}_{\perp }\equiv \mathcal {P}_\perp \boldsymbol {f}$ with $\boldsymbol {f}=[f_\mathcal {r},f_\theta,f_\varphi ]^{{\rm T}}$ and projection matrix $\mathcal {P}_\perp \equiv ( \begin {smallmatrix}0 & 1 & 0 \\ 0 & 0 & 1 \end {smallmatrix})$, and the double square brackets denote the difference between the top and bottom values evaluated on the inner sides of the bounding surfaces,

(2.8)\begin{equation} [\![ {\phi} ]\!]\equiv\phi|_{\eta_{T^-}}-\phi|_{\eta_{B^+}}, \quad \text{with } \phi|_{\eta_{T^-}}\equiv\lim_{\mathcal{r}\to\eta_{T^-}} \phi \quad \text{and} \quad \phi|_{\eta_{B^+}}\equiv\lim_{\mathcal{r}\to\eta_{B^+}} \phi, \end{equation}

and the in-plane divergence operator is

(2.9)\begin{equation} {\nabla}_{{\perp}}^\mathcal{R}\boldsymbol{\cdot}(\phi{\boldsymbol{f}}_{{\perp}}) =\frac{1}{\mathcal{R}\sin\theta}\left(\frac{\partial(\phi f_\theta\sin\theta)}{\partial\theta}+\frac{\partial(\phi f_\varphi)}{\partial\varphi}\right). \end{equation}

Using (2.3), (2.6) and (2.7), the depth average of the governing equations (2.1) is

(2.10)\begin{align} &\frac{\partial}{\partial t} \left[\begin{matrix} \mathcal{H} \langle \rho \rangle \\ \mathcal{H} \langle \rho v_\mathcal{r} \rangle \\ \mathcal{H} \langle \rho {\boldsymbol{v}}_{{\perp}} \rangle \\ \mathcal{H} \langle \rho e_{T} \rangle \end{matrix}\right] + {\nabla}_{{\perp}}^\mathcal{R} \boldsymbol{\cdot} \left[\begin{matrix} \mathcal{L}\mathcal{R} \overline{\rho{\boldsymbol{v}}_{{\perp}} } \\ \mathcal{L}\mathcal{R} \overline{\rho v_\mathcal{r}\,{\boldsymbol{v}}_{{\perp}} } \\ \mathcal{L}\mathcal{R}(\overline{\rho{\boldsymbol{v}}_{{\perp}} \otimes{\boldsymbol{v}}_{{\perp}} } + \bar{p}{\boldsymbol{\mathsf{I}}}_2) \\ \mathcal{L}\mathcal{R}(\overline{\rho e_{T}{\boldsymbol{v}}_{{\perp}} } + \overline{p{\boldsymbol{v}}_{{\perp}} }) \end{matrix}\right]\nonumber\\ &\quad = \left[\begin{matrix} 0 \\ -[\![ {p} ]\!]- \langle \rho \rangle g \mathcal{H} \\ [\![ {p{\nabla}_{{\perp}}^\eta\mathscr{Z}} ]\!] \\ [\![ {p\boldsymbol{v}\boldsymbol{\cdot}\nabla^\eta\mathscr{Z}} ]\!]- \langle \rho v_\mathcal{r} \rangle g\mathcal{H} \end{matrix}\right] -2\mathcal{L} \left[\begin{matrix} \overline{\rho v_\mathcal{r}} \\ \overline{\rho v_\mathcal{r} v_\mathcal{r}} \\ \overline{\rho v_\mathcal{r}{\boldsymbol{v}}_{{\perp}} } \\ \overline{\rho v_\mathcal{r} e_T}+\overline{p v_\mathcal{r}} \end{matrix}\right]{,} \end{align}

with $g \equiv \boldsymbol {g}\boldsymbol {\cdot }\boldsymbol {e}_\mathcal {r}$ (note that $\boldsymbol {g}\boldsymbol {\cdot }\boldsymbol {e}_\theta = \boldsymbol {g}\boldsymbol {\cdot }\boldsymbol {e}_\varphi = 0$ from the assumptions).

2.3. Thin-layer assumption

Let $d$ denote the characteristic length of the layer depth. Assuming the layer is thin compared with its radius ($d\ll \mathcal {R}$), with $\mathcal {R}\equiv \mathcal {R}^*/\ell _o$, where $\mathcal {R}^*$ is the characteristic radius of the planet (for Earth $\mathcal {R}^*=6371\ \mbox {km}$). Defining $\delta \equiv d/\mathcal {R}$ and using the results of Appendix B.2 ($\bar {\phi }=\langle \phi \rangle +\mathcal {O}(\delta ^2)$), (2.6) and (2.7) become

(2.11)\begin{gather} \int_{\eta_B}^{\eta_T}\frac{\partial \phi}{\partial t}\, \mathrm{d} \mathcal{r}=\frac{\partial }{\partial t}(\mathcal{H} \bar{\phi}) - \left[\!\!\left[ {\phi\frac{\partial\eta}{\partial t}} \right]\!\!\right]+\mathcal{O}(\delta^2), \end{gather}

(2.12)\begin{gather}\int_{\eta_B}^{\eta_T}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{f} \, \mathrm{d} \mathcal{r}= {\nabla}_{{\perp}}^\mathcal{R}\boldsymbol{\cdot}(\mathcal{H}{\bar{\boldsymbol{f}}}_{{\perp}}) +[\![ \,{\boldsymbol{f}\boldsymbol{\cdot}\boldsymbol{n}} ]\!]+ 2\frac{\mathcal{H}}{\eta_B}(\,\bar{\boldsymbol{f}}\boldsymbol{\cdot}\boldsymbol{e}_\mathcal{r}+ [\![ \, {{\boldsymbol{f}}_{{\perp}}\boldsymbol{\cdot}{\boldsymbol{n}}_{{\perp}}} ]\!]) +\mathcal{O}(\delta^2), \end{gather}

with the normal vector defined as $\boldsymbol {n}_i\equiv -\nabla ^{\mathcal {R}}\mathscr {Z}_i$ and ${\boldsymbol {n}_i}_{\perp }={\mathcal {P}}_{\perp }\boldsymbol {n}_i$ (see Appendix B.2 for details).

Applying the above and using the surface evolution equation (2.3) on (2.10) results in the depth average thin-layer equations

(2.13)\begin{align} &\frac{\partial}{\partial t} \left[\begin{matrix} \mathcal{H} \bar{\rho} \\ \mathcal{H} \overline{\rho v_\mathcal{r}} \\ \mathcal{H} \overline{\rho {\boldsymbol{v}}_{{\perp}} } \\ \mathcal{H} \overline{\rho e_{T}} \end{matrix}\right] + {\nabla}_{{\perp}}^\mathcal{R} \boldsymbol {\cdot} \left[\begin{matrix} \mathcal{H} \overline{\rho{\boldsymbol{v}}_{{\perp}} } \\ \mathcal{H} \overline{\rho v_\mathcal{r}\,{\boldsymbol{v}}_{{\perp}} } \\ \mathcal{H} \overline{\rho{\boldsymbol{v}}_{{\perp}} \otimes{\boldsymbol{v}}_{{\perp}} }+ \mathcal{H} \bar{p}{\boldsymbol{\mathsf{I}}}_2 \\ \mathcal{H} \overline{\rho e_{T}{\boldsymbol{v}}_{{\perp}} }+ \mathcal{H} \overline{p{\boldsymbol{v}}_{{\perp}} } \end{matrix}\right]\nonumber\\ &\quad = \left[\begin{matrix} 0 \\ -[\![ {p} ]\!]- \bar{\rho}g \mathcal{H} \\ -[\![ {p{\boldsymbol{n}}_{{\perp}}} ]\!] \\ -[\![ {p\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n}} ]\!]- \overline{\rho v_\mathcal{r}}g\mathcal{H} \end{matrix}\right] -2\frac{\mathcal{H}}{\eta_B} \left[\begin{matrix} \overline{\rho v_\mathcal{r}} \\ \overline{\rho v_\mathcal{r} v_\mathcal{r}} \\ \overline{\rho v_\mathcal{r}{\boldsymbol{v}}_{{\perp}} } + [\![ {p{\boldsymbol{n}}_{{\perp}}} ]\!]/2 \\ \overline{\rho v_\mathcal{r} e_T}+\overline{p v_\mathcal{r}}+[\![ {p{\boldsymbol{v}}_{{\perp}} \boldsymbol{\cdot}{\boldsymbol{n}}_{{\perp}}} ]\!]/2 \end{matrix}\right] + \mathcal{O}(\delta^2). \end{align}

2.4. Long-wave assumption

Considering a perturbation to the top surface $\eta _T$ with a wavelength $L$ and amplitude $a$, the relative water depth is defined as $\varepsilon \equiv d/L$, where $d$ is the characteristic vertical length and $a\ll d$. Then, assuming long waves or shallow water ($\varepsilon \ll 1$) the vertical velocity and the gradient of the surface become $\bar {v}_\mathcal {r}=\mathcal {O}(\varepsilon )$ and ${\nabla }_{\perp }^\mathcal {R}\eta =\mathcal {O}(a/L)$, and (2.13) becomes

(2.14)\begin{align} &\frac{\partial}{\partial t} \left[\begin{matrix} \mathcal{H} \bar{\rho} \\ \mathcal{H} \overline{\rho v_\mathcal{r}} \\ \mathcal{H} \overline{\rho {\boldsymbol{v}}_{{\perp}} } \\ \mathcal{H} \overline{\rho e_{T}} \end{matrix}\right] + {\nabla}_{{\perp}}^\mathcal{R} \boldsymbol{\cdot} \left[\begin{matrix} \mathcal{H} \overline{\rho{\boldsymbol{v}}_{{\perp}} } \\ \mathcal{H} \overline{\rho v_\mathcal{r} {\boldsymbol{v}}_{{\perp}} } \\ \mathcal{H} \overline{\rho{\boldsymbol{v}}_{{\perp}} \otimes{\boldsymbol{v}}_{{\perp}} } + \mathcal{H} \bar{p}{\boldsymbol{\mathsf{I}}}_2 \\ \mathcal{H} \overline{\rho e_{T}{\boldsymbol{v}}_{{\perp}} } + \mathcal{H} \overline{p{\boldsymbol{v}}_{{\perp}} } \end{matrix}\right]\nonumber\\ &\quad = \left[\begin{matrix} 0 \\ -[\![ {p} ]\!]- \bar{\rho}g\mathcal{H} \\ -[\![ {p{\boldsymbol{n}}_{{\perp}}} ]\!] \\ -[\![ {p\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n}} ]\!] - \overline{\rho v_\mathcal{r}}g\mathcal{H} \end{matrix}\right] + \mathcal{O}(\delta\varepsilon), \end{align}

and the surface evolution equation (2.3) becomes

(2.15)\begin{equation} \frac{\partial\eta_i}{\partial t}-\boldsymbol{v}|_{\eta_i}\boldsymbol{\cdot}\boldsymbol{n}_i=\mathcal{O}(\delta\varepsilon).\end{equation}

2.5. Density-weighted average decomposition

Similarly to common practice in compressible-turbulence studies, density-weighted averages are recast following a Favre decomposition, $\phi =\tilde {\phi }+\phi ''$, where $\tilde {\phi }=\overline {\rho \phi }/\bar {\rho }$. Here, the Reynolds decomposition uses the logarithmic depth average as $\phi =\bar {\phi }+\phi '$. With these decompositions, (2.14) becomes

(2.16)\begin{gather} \frac{\partial}{\partial t} \left[ \begin{matrix} \mathcal{H} \bar{\rho} \\ \mathcal{H} \bar{\rho}\tilde{v}_\mathcal{r} \\ \mathcal{H} \bar{\rho}{\tilde{\boldsymbol{v}}}_{{\perp}} \\ \mathcal{H} \bar{\rho}\tilde{e}_T \end{matrix} \right] \!+\! {\nabla}_{{\perp}}^\mathcal{R} \boldsymbol{\cdot} \left[ \begin{matrix} \mathcal{H} \bar{\rho}{\tilde{\boldsymbol{v}}}_{{\perp}} \\ \mathcal{H} \bar{\rho}\tilde{v}_\mathcal{r}{\tilde{\boldsymbol{v}}}_{{\perp}} \\ \mathcal{H} \bar{\rho}{\tilde{\boldsymbol{v}}}_{{\perp}}\otimes{\tilde{\boldsymbol{v}}}_{{\perp}} \!+\! \mathcal{H}\bar{p}{\boldsymbol{\mathsf{I}}}_2 \\ \mathcal{H} \bar{\rho}\tilde{e}_T{\tilde{\boldsymbol{v}}}_{{\perp}} + \mathcal{H} \bar{p}{\tilde{\boldsymbol{v}}}_{{\perp}} \end{matrix} \right] = \left[ \begin{matrix} 0 \\ -[\![ {p} ]\!]- \bar{\rho}g\mathcal{H} \\ -[\![ {p{\boldsymbol{n}}_{{\perp}}} ]\!] \\ -[\![ {p\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n}} ]\!] - \bar{\rho}\tilde{v}_{\mathcal{r}} g\mathcal{H} \end{matrix} \right] + \mathcal{C} + \mathcal{O}(\delta\varepsilon),\end{gather}

where the in-plane velocity is decomposed as ${\boldsymbol {v}}_{\perp } ={\tilde {\boldsymbol {v}}}_{\perp }+{\boldsymbol {v}}_{\perp }''$. The closure vector $\mathcal {C}$ corresponds to

(2.17)\begin{align} \mathcal{C}&={-}\partial \left[ \begin{matrix} 0, & 0, & 0, & \mathcal{H}\bar{\rho}\widetilde{{\boldsymbol{v}}_{{\perp}}^{''} \boldsymbol{\cdot}{\boldsymbol{v}}_{{\perp}}^{''}}/2 \end{matrix} \right]^{{\rm T}}/\partial t \nonumber\\ & \quad - {\nabla}_{{\perp}}^\mathcal{R} \boldsymbol {\cdot} \left[ \begin{matrix} 0, & \mathcal{H} \bar{\rho}\widetilde{v_\mathcal{r}^{\prime\prime}{\boldsymbol{v}}_{{\perp}}^{\prime\prime}}, & \mathcal{H} \bar{\rho}\widetilde{{\boldsymbol{v}}_{{\perp}}^{\prime\prime}\otimes{\boldsymbol{v}}_{{\perp}}^{\prime\prime}}, & \mathcal{H} \bar{\rho}\widetilde{ e_T^{\prime\prime}{\boldsymbol{v}}_{{\perp}}^{\prime\prime}}+ \mathcal{H} (\overline{p'{\boldsymbol{v}'}_{{\perp}}}+\bar{p}\overline{{\boldsymbol{v}}_{{\perp}}^{\prime\prime}}) \end{matrix} \right]^{{\rm T}}. \end{align}

Using Taylor series expansions, the closure vector $\mathcal {C}$ is proven to be $\mathcal {O}(\delta \varepsilon )$ or, as shown in (B23), higher (see Appendix B.3 for details) and (2.16) becomes

(2.18)\begin{equation} \frac{\partial}{\partial t} \left[ \begin{matrix} \mathcal{H} \bar{\rho} \\ \mathcal{H} \bar{\rho}\tilde{v}_\mathcal{r} \\ \mathcal{H} \bar{\rho}{\tilde{\boldsymbol{v}}}_{{\perp}} \\ \mathcal{H} \bar{\rho}\tilde{e}_T \end{matrix} \right] + {\nabla}_{{\perp}}^\mathcal{R} \boldsymbol {\cdot} \left[ \begin{matrix} \mathcal{H} \bar{\rho}{\tilde{\boldsymbol{v}}}_{{\perp}} \\ \mathcal{H} \bar{\rho}\tilde{v}_\mathcal{r}{\tilde{\boldsymbol{v}}}_{{\perp}} \\ \mathcal{H} \bar{\rho}{\tilde{\boldsymbol{v}}}_{{\perp}}\otimes{\tilde{\boldsymbol{v}}}_{{\perp}} + \mathcal{H}\bar{p}{\boldsymbol{\mathsf{I}}}_2 \\ \mathcal{H} \bar{\rho}\tilde{e}_T{\tilde{\boldsymbol{v}}}_{{\perp}} + \mathcal{H} \bar{p}{\tilde{\boldsymbol{v}}}_{{\perp}} \end{matrix} \right] = \left[ \begin{matrix} 0 \\ -[\![ {p} ]\!]- \bar{\rho}g\mathcal{H} \\ -[\![ {p{\boldsymbol{n}}_{{\perp}}} ]\!] \\ -[\![ {p\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n}} ]\!] - \bar{\rho}\tilde{v}_{\mathcal{r}} g\mathcal{H} \end{matrix} \right] + \mathcal{O}(\delta\varepsilon). \end{equation}

2.6. Long-wave expansion

The independent variables are rescaled based on the relative water depth $\varepsilon$ as

(2.19)\begin{equation} t=\varepsilon^{{-}1}t^o,\quad {\boldsymbol{x}}_{{\perp}}=\varepsilon^{{-}1}{\boldsymbol{x}}_{{\perp}}^o,\quad \mathcal{r}=\mathcal{r}^o,\end{equation}

where the superscript $^o$ denotes the rescaled quantities, and ${\boldsymbol {x}}_{\perp }$ is a position vector on the plane formed by the $\boldsymbol {e}_\theta$ and $\boldsymbol {e}_\varphi$ vectors. With this rescaling the variables become

(2.20)\begin{equation} \mathcal{H} = \mathcal{H}^o ,\quad \bar{\rho} = \bar{\rho}^o ,\quad \bar{p} = \bar{p}^o ,\quad \tilde{v}_\mathcal{r} = \varepsilon \tilde{v}_\mathcal{r}^o ,\quad {\tilde{\boldsymbol{v}}}_{{\perp}} = {\tilde{\boldsymbol{v}}}_{{\perp}}^o ,\quad \tilde{e}_T = \tilde{e}_T^o . \end{equation}

Analogous to standard practices in the derivation of the SWE (see, e.g. Narayanan Reference Narayanan2003), the coefficients of the rescaled variables are based on a set of modelling choices over the resulting equations. For the purposes of this work, these coefficients are derived based on the following constraints over the resulting leading-order equations and rescaled variables.

1. (i) The vertical momentum equation corresponds to the hydrostatic balance.

2. (ii) All the terms in the continuity equation are of the same order.

3. (iii) All the terms in the in-plane momentum equation are of the same order.

4. (iv) Pressure and density changes follow the same scaling (i.e. they are related via the isentropic constraint).

Constraints (i)–(iii) are used to recover the SWE in the incompressible-flow case and constraint (iv) is added for compatibility with the internal-energy equation in the compressible-flow case. Note that the resulting leading-order equations will differ for another set of considerations. Appendix B.4 describes the derivation of the coefficients and the leading-order equations. From here on, the superscript $^o$ is omitted for readability, and the resulting leading-order equations are

(2.21)\begin{align} \frac{\partial}{\partial t} \left[ \begin{matrix} \mathcal{H} \bar{\rho} \\ 0 \\ \mathcal{H} \bar{\rho}{\tilde{\boldsymbol{v}}}_{{\perp}} \\ \mathcal{H} \bar{\rho}\tilde{e}_T \end{matrix} \right] +{\nabla}_{{\perp}}^{\mathcal{R}} \boldsymbol{\cdot} \left[ \begin{matrix} \mathcal{H} \bar{\rho}{\tilde{\boldsymbol{v}}}_{{\perp}} \\ 0 \\ \mathcal{H} \bar{\rho}{\tilde{\boldsymbol{v}}}_{{\perp}}\otimes{\tilde{\boldsymbol{v}}}_{{\perp}} + \mathcal{H}\bar{p}{\boldsymbol{\mathsf{I}}}_2 \\ \mathcal{H} \bar{\rho}\tilde{e}_T{\tilde{\boldsymbol{v}}}_{{\perp}} + \mathcal{H} \bar{p}{\tilde{\boldsymbol{v}}}_{{\perp}} \end{matrix} \right] = \left[ \begin{matrix} 0 \\ -[\![ {p} ]\!]- \bar{\rho}g\mathcal{H} \\ -[\![ {p{\boldsymbol{n}}_{{\perp}}} ]\!] \\ -[\![ {p\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n}} ]\!] - \bar{\rho}\tilde{v}_{\mathcal{r}} g\mathcal{H} \end{matrix} \right] + \mathcal{O}(\delta\varepsilon), \end{align}

and

(2.22)\begin{equation} \frac{\partial\eta_i}{\partial t}-\boldsymbol{v}|_{\eta_i}\boldsymbol{\cdot}\boldsymbol{n}_i=\mathcal{O}(\delta\varepsilon). \end{equation}

2.7. Isentropic-flow constraint

The kinetic energy equation is obtained by taking the dot product of the Favre average velocity $\tilde {\boldsymbol {v}}$ with the momentum-balance statement in (2.21), including the radial components. Subtracting the result from the total energy equation in (2.21) yields the governing equation for the Favre-averaged internal energy after application of the continuity equation ${\nabla }_{\perp }^\mathcal {R}\boldsymbol {\cdot }{\tilde {\boldsymbol {v}}}_{\perp }= -(\mathcal {H}\bar {\rho })^{-1}{\rm D}(\mathcal {H}\bar {\rho })/{\rm D}t$, i.e.

(2.23)\begin{equation} \mathcal{H}\bar{\rho}\frac{{\rm D} \tilde{e}}{{\rm D} t} -\frac{\bar{p}}{\bar{\rho}}\frac{{\rm D} (\mathcal{H}\bar{\rho})}{{\rm D} t} +[\![ {p\boldsymbol{v}''\boldsymbol{\cdot}\boldsymbol{n}} ]\!]=\mathcal{O}(\delta\varepsilon).\end{equation}

Rearranging the material derivative of $\mathcal {H}\bar {\rho }$, introducing $p = \bar {p} + p'$, noting that $[\![ {\bar {p}\boldsymbol {f}} ]\!] = \bar {p}[\![ {\boldsymbol {f}} ]\!]$ and $[\![ {p'\boldsymbol {v}''\boldsymbol {\cdot }\boldsymbol {n}} ]\!]=\mathcal {O}(\delta \varepsilon )$ (see Appendix B.3 for details) gives

(2.24)\begin{equation} \mathcal{H}\bar{\rho}\frac{{\rm D} \tilde{e}}{{\rm D}t}={-}\mathcal{H}\bar{\rho}\,\bar{p}\frac{{\rm D}}{{\rm D}t}\left(\frac{1}{\bar{\rho}}\right) + \bar{p}\left( \frac{{\rm D} \mathcal{H}}{{\rm D}t} - [\![ {\boldsymbol{v}''\boldsymbol{\cdot}\boldsymbol{n}} ]\!]\right) + \mathcal{O}(\delta\varepsilon). \end{equation}

Applying (2.22) at $\eta _T$ and $\eta _B$, introducing the Favre decomposition ${\boldsymbol {v}}_{\perp } = {\tilde {\boldsymbol {v}}}_{\perp } + {\boldsymbol {v}}_{\perp }''$ and $v_\mathcal {r} = \tilde {v}_\mathcal {r} + v_\mathcal {r}''$, and subtracting the results yields $D\mathcal {H}/Dt = [\![ {v_\mathcal {r}''} ]\!] + \mathcal {O}(\delta \varepsilon )$. The internal-energy equation becomes

(2.25)\begin{equation} \mathcal{H}\bar{\rho}\frac{{\rm D} \tilde{e}}{{\rm D}t}={-}\mathcal{H}\bar{\rho}\,\bar{p}\frac{{\rm D}}{{\rm D}t}\left(\frac{1}{\bar{\rho}}\right) + \bar{p}[\![ {\boldsymbol{v}''\boldsymbol{\cdot}(\boldsymbol{e}_\mathcal{r}-\boldsymbol{n})} ]\!] + \mathcal{O}(\delta\varepsilon). \end{equation}

The long-wave assumption implies that $\boldsymbol {n} = \boldsymbol {e}_\mathcal {r} + \mathcal {O}(\varepsilon )$. The thin-layer assumption implies that $\boldsymbol {v}'' = \mathcal {O}(\delta )$. Hence, $\boldsymbol {v}''\boldsymbol {\cdot }(\boldsymbol {e}_\mathcal {r}-\boldsymbol {n}) = \mathcal {O}(\delta \varepsilon$), and the internal-energy equation reduces to

(2.26)\begin{equation} \mathcal{H}\bar{\rho}\left[\frac{{\rm D} \tilde{e}}{{\rm D}t}+\bar{p}\frac{{\rm D}}{{\rm D}t}\left(\frac{1}{\,\bar{\rho}\,}\right)\right] = \mathcal{O}(\delta\varepsilon). \end{equation}

Equation (2.26) implies that, to first order in $\delta \varepsilon$, the entropy of a material element of a thin layer is conserved along its trajectory in space–time. The flow is therefore isentropic.

Let $\mathcal {V}\equiv 1/\bar {\rho }$ be the specific volume of the layer (not necessarily the depth average value of $\rho ^{-1}$). Let $\mathcal {S}$ be the specific entropy of the layer (not necessarily the Favre-averaged value of the specific entropy). Considering that $\tilde {e}$ depends on the state variables $\mathcal {V}$ and $\mathcal {S}$ (having previously ignored any chemical-potential dependency), the total differential of $\tilde {e}$ is

(2.27)\begin{equation} \mathrm{d}\tilde{e}=\left(\frac{\partial\tilde{e}}{\partial \mathcal{V}}\right)_{\mathcal{S}}\, \mathrm{d}\mathcal{V}+\left(\frac{\partial\tilde{e}}{\partial\mathcal{S}} \right)_{\mathcal{V}}\, \mathrm{d}\mathcal{S}. \end{equation}

Assuming that any material element of the layer evolves under the isentropic constraint $\mathrm {d}\mathcal {S}=0$, (2.26) and (2.27) give

(2.28)\begin{equation} \bar{p} ={-} \left(\frac{\partial\tilde{e}}{\partial \mathcal{V}}\right)_{\mathcal{S}}. \end{equation}

Since the thermodynamic pressure $\bar {p}$ depends on the state variables $\mathcal {V}$ and $\mathcal {S}$, the isentropic constraint imposes that

(2.29)\begin{equation} \mathrm{d}\bar{p}=\left(\frac{\partial\bar{p}}{\partial \mathcal{V}} \right)_{\mathcal{S}}\, \mathrm{d}\mathcal{V} = \left(\frac{\partial\bar{p}}{\partial \bar{\rho}}\right)_{\mathcal{S}}\, \mathrm{d}\bar{\rho} = \mathscr{C}^2\, \mathrm{d}\bar{\rho}, \end{equation}

where $\mathscr {C} \equiv [(\partial \bar {p}/\partial \bar {\rho })_{\mathcal {S}}]^{1/2}$ is the isentropic sound speed of the layer.

The total energy equation in (2.21) can thus be replaced by

(2.30)\begin{equation} \frac{{\rm D} \bar{p}}{{\rm D} t} = \mathscr{C}^2 \frac{{\rm D} \bar{\rho}}{{\rm D} t}, \end{equation}

where $\mathscr {C}$ is evaluated from the depth-averaged variables.

3.1. Shallow ocean equations

The seawater density, $\rho _w$, is assumed constant and uniform. The continuity equation in (2.21) simplifies to ${\rm D}^{\boldsymbol {U}} H/{\rm D} t = - H {\nabla }_{\perp }^\mathcal {R}\boldsymbol {\cdot }\boldsymbol {U}$. The radial component of the momentum equation in (2.21) becomes $[\![ {p} ]\!] = p|_{\eta _1} - p|_{\eta _0} = -\rho _wg H$. This is the hydrostatic balance (‘$\mathrm {d} p = - \rho _wg\,\mathrm {d} \mathcal {r}$’), which, in the uniform-density case integrates to the linear profile $p(\mathcal {r}) = p|_{\eta _1} - \rho _wg(\mathcal {r}-\eta _1)$. Depth averaging gives $\bar {p}=p|_{\eta _1} + \rho _wg H/2 +\mathcal {O}(\delta ^2)$. Substituting these in (2.21) yields

(3.4)\begin{equation} \frac{{\rm D}^{\boldsymbol{U}}}{{\rm D} t} \left[\begin{matrix} H \\ \boldsymbol{U} \end{matrix}\right] ={-}\left[\begin{matrix} H {\nabla}_{{\perp}}^\mathcal{R}\boldsymbol{\cdot}\boldsymbol{U} \\ \rho_w^{{-}1}{\nabla}_{{\perp}}^\mathcal{R} p|_{\eta_1} + g{\nabla}_{{\perp}}^\mathcal{R}\eta_1 \end{matrix}\right].\end{equation}

Note that the internal-energy equation is removed from the ocean-layer dynamics since ‘pressure’ in the uniform-density (and isothermal) framework purely assumes a kinematic role and not a thermodynamic one. Retaining the internal-energy equation would over-constrain the dynamical system.

Equation (3.4) is strictly the same dynamical equation as that of the OWC models (see (1.1)). However, note that the TWC model will let $\eta _1$, $\boldsymbol {U}$ and $p|_{\eta _1}$ all be influenced by the atmospheric layer, whereas OWC models prescribe assumed space–time dependencies for $p|_{\eta _1}$ that do not depend on the evolution of the system. In the absence of a pressure gradient on the ocean surface and a fixed seabed position, (3.4) is in a closed form and is strictly equivalent to the classical SWE (see, e.g. Vallis Reference Vallis2017), referred here to as the zero-way coupling (ZWC) model.

3.2. Shallow atmosphere equations

The continuity and in-plane momentum equations in (2.21) are rearranged in primitive form and the total energy equation replaced by (2.30) to give

(3.5)\begin{equation} \frac{{\rm D} ^{\boldsymbol{u}}}{{\rm D} t} \left[\begin{matrix} \varrho \\ \boldsymbol{u} \\ {\rm \pi}\end{matrix}\right] ={-}\left[\begin{matrix} \varrho\varPsi \\ (\varrho h)^{{-}1}({\nabla}_{{\perp}}^\mathcal{R}(h{\rm \pi})-p|_{\eta_1} {\nabla}_{{\perp}}^\mathcal{R} h) -g{\nabla}_{{\perp}}^\mathcal{R}\eta_2 \\ \varrho\varPsi\mathscr{C}^2 \end{matrix}\right],\end{equation}

where $\varPsi \equiv h^{-1}[{\nabla }_{\perp }^\mathcal {R}\boldsymbol {\cdot }(h\boldsymbol {u})+\partial h/\partial t]$ and $p|_{\eta _1} = p|_{\eta _2} + \varrho g h$.

Air is assumed to be thermally ($p^* = \rho ^* \mathscr {R}_{air} T^*$) and calorically ($\mathrm {d}e^* = c_v \, \mathrm {d}T^*$) perfect (humidity and phase changes are not taken into account), where the gas constant $\mathscr {R}_{air}$ and the specific heat at constant volume $c_v$ are true dimensional constants. Further assuming $e^*(T^*=0)=0$ results in $e^*=c_v T^*$. Using the reference pressure $p_o = \rho _o u_o^2$ and defining the non-dimensional temperature as $T \equiv \mathscr {R}_{air} T^*/u_o^2$, the thermal equation of state (EoS) becomes $p=\rho T$ and the caloric EoS becomes $e = (c_v/\mathscr {R}_{air}) T$. Depth-averaging both gives $\bar {p} = \bar {\rho }\tilde {T}$ and $\tilde {e} = (c_v/\mathscr {R}_{air}) \tilde {T}$, and injecting these in (2.27) using (2.28) yields

(3.6)\begin{equation} \mathrm{d}\bar{p}=\gamma\frac{\bar{p}}{\bar{\rho}}\, \mathrm{d}\bar{\rho} +\bar{\rho}\frac{\mathscr{R}_{air}}{c_v} \left(\frac{\partial\tilde{e}}{\partial\mathcal{S}}\right)_{\mathcal{V}}\, \mathrm{d}\mathcal{S}, \end{equation}

where $\gamma \equiv c_p/c_v$, $c_p$ is the specific heat at constant pressure, and $\mathscr {R}_{air} = c_p-c_v$ (Mayer's relation). Thus, using (2.29), the isentropic sound speed of the layer is

(3.7)\begin{equation} \mathscr{C} = \sqrt{\gamma \frac{\bar{p}}{\bar{\rho}}} = \sqrt{\gamma \frac{\rm \pi}{\varrho}}, \end{equation}

where ${\rm \pi}$ and $\varrho$ are, respectively, the (logarithmic) depth-averaged pressure and density. Section 5 will show that in standard atmospheric conditions this speed, evaluated on the entire thickness of the troposphere (and beyond), matches that of the so-called Lamb wave invoked in OWC models.

3.3. Boundary conditions

3.3.1. Edge boundary conditions

The seabed ($\eta _0$) is non-uniform but stationary. The outer-space boundary ($\eta _2$) is assumed uniform and stationary. The steady edge boundary conditions therefore impose that $H + h = \eta _2-\eta _0$ remains constant with time,

(3.8)\begin{equation} \frac{\partial H}{\partial t}+\frac{\partial h}{\partial t}=0.\end{equation}

The rather strict constraints on $\eta _2$ can be relaxed in future studies by adding further atmospheric layer(s) to better account for the high-atmosphere physics (Aa et al. Reference Aa, Zhang, Erickson, Vierinen, Coster, Goncharenko, Spicher and Rideout2022), which is beyond the scope of this study, where $\eta _2$ is placed at the interface with the ionosphere, which is modelled as a uniform semi-infinite vacuum (i.e. Earth is placed in a vacuum from the ionosphere outwards).

3.3.2. Interface conditions

The surface evolution equation (2.22) is evaluated on the upper ($+$) and lower ($-$) sides of the $i^\mathrm {th}$ layer,

(3.9)\begin{equation} \frac{\partial\eta_i^{{\pm}}}{\partial t} = \boldsymbol{v}|_{\eta_i^{{\pm}}}\boldsymbol{\cdot}\boldsymbol{n}_i^{{\pm}}+\mathcal{O}(\delta\varepsilon). \end{equation}

It follows that two neighbouring layers remain in contact if $\eta _i^+ = \eta _i^{-}$ at all times, i.e. if $\boldsymbol {v}|_{\eta _i^+}\boldsymbol {\cdot }\boldsymbol {n}_i^+ = \boldsymbol {v}|_{\eta _i^-}\boldsymbol {\cdot }\boldsymbol {n}_i^-$. Since $\boldsymbol {v}$ (the local velocity field) is continuous, the surface evolution equation enforces that no layer detaches from its neighbour. The ocean and atmosphere are thus in contact at all times by simply requiring that they initially satisfy  $\eta _1^+(t=0) = \eta _1^-(t=0) = \eta _1(t=0)$. Note, however, that depth-averaged velocities are not required (nor expected) to be continuous across the interface (e.g. the depth-averaged velocity of the atmosphere will generally not be equal to that of the ocean).

3.3.3. Pressure boundary conditions

Pressure is considered continuous across the ocean-atmosphere interface (e.g. surface tension is negligible on long waves). On the outer-space side, pressure is assumed to be zero given the vacuum assumption i.e. $p|_{\eta _2}=0$. The jump condition (from the radial component of the governing equation) across the atmospheric layer simplifies to

(3.10)\begin{equation} p|_{\eta_1} = \varrho g h.\end{equation}

3.4. Two-way coupled equations

Combining (3.8) with the water-layer continuity equation (3.4) yields

(3.11)\begin{equation} \frac{\partial h}{\partial t}= {\nabla}_{{\perp}}^\mathcal{R}\boldsymbol{\cdot}(H \boldsymbol{U}). \end{equation}

Moreover, boundary conditions $p|_{\eta _2} = 0$ and ${\nabla }_{\perp }^\mathcal {R}\eta _2 = \boldsymbol {0}$ reduce the atmosphere-layer equation to

(3.12)\begin{equation} \frac{{\rm D} ^{\boldsymbol{u}}}{{\rm D} t} \left[\begin{matrix} \varrho \\ \boldsymbol{u} \\ {\rm \pi}\end{matrix}\right] ={-}\left[\begin{matrix} \varrho\varPsi \\ (\varrho h)^{{-}1}{\nabla}_{{\perp}}^\mathcal{R}(h{\rm \pi}) +g{\nabla}_{{\perp}}^\mathcal{R}\eta_1\\ \varrho\varPsi\mathscr{C}^2 \end{matrix}\right],\end{equation}

with $\varPsi = h^{-1}{\nabla }_{\perp }^\mathcal {R}\boldsymbol {\cdot }(h\boldsymbol {u}+H\boldsymbol {U})$.

Equations (3.4) and (3.12) can be solved simultaneously to form the TWC model

(3.13) \begin{equation} \frac{\partial}{\partial t} \left[\begin{matrix} \eta_1 \\ \boldsymbol{U} \\ \varrho \\ \boldsymbol{u} \\ {\rm \pi} \end{matrix}\right] + \left[\begin{matrix} \boldsymbol{U}\boldsymbol{\cdot}{\nabla}_{{\perp}}^\mathcal{R} H \\ (\boldsymbol{U}\boldsymbol{\cdot}{\nabla}_{{\perp}}^\mathcal{R})\boldsymbol{U} \\ \boldsymbol{u}\boldsymbol{\cdot}{\nabla}_{{\perp}}^\mathcal{R} \varrho \\ (\boldsymbol{u}\boldsymbol{\cdot}{\nabla}_{{\perp}}^\mathcal{R})\boldsymbol{u} \\ \boldsymbol{u}\boldsymbol{\cdot}{\nabla}_{{\perp}}^\mathcal{R} {\rm \pi} \end{matrix}\right] ={-}\left[\begin{matrix} H{\nabla}_{{\perp}}^\mathcal{R}\boldsymbol{\cdot}\boldsymbol{U} \\ \rho_w^{{-}1}{\nabla}_{{\perp}}^\mathcal{R}(\varrho g h) + g{\nabla}_{{\perp}}^\mathcal{R}\eta_1 \\ \varrho\varPsi \\ (\varrho h)^{{-}1}{\nabla}_{{\perp}}^\mathcal{R}(h{\rm \pi}) +g{\nabla}_{{\perp}}^\mathcal{R}\eta_1 \\ \varrho\varPsi\mathscr{C}^2 \end{matrix}\right],\end{equation}

with $h = \eta _2-\eta _1$, $H = \eta _1-\eta _0$, $\varPsi = h^{-1}{\nabla }_{\perp }^\mathcal {R}\boldsymbol {\cdot }(h\boldsymbol {u}+H\boldsymbol {U})$, $\mathscr {C}^2 = \gamma {\rm \pi}/\varrho$, where $\gamma$ is a constant (set to $1.4$). Note that $\eta _0$ and $\eta _2$ are set initially and do not change, hence, the choice to integrate $\eta _1$ rather than $H$ and $h$ that can be evaluated directly from $\eta _1$ and the initial $\eta _0$ and $\eta _2$ profiles.

4.1. Quasi-linear one-dimensional form

The one-dimensional quasi-linear form is found by projecting (3.13) along a great circle with coordinate $s$ along which $\boldsymbol {e}_{s}$ is the unit vector tangent to the sphere of radius $\mathcal {R}$. Letting $\boldsymbol {\mu } \equiv [\eta _1,U,\varrho,u,{\rm \pi} ]^{{\rm T}}$ denote the vector of primitive variables, with $U=\boldsymbol {U}\boldsymbol {\cdot }\boldsymbol {e}_{s}$ and $u=\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {e}_{s}$, the projected equation is

(4.1)\begin{equation} \frac{\partial\boldsymbol{\mu}}{\partial t} + {\boldsymbol{\mathsf{A}}} \frac{\partial\boldsymbol{\mu}}{\partial s} =\boldsymbol{0}, \quad\text{with } {\boldsymbol{\mathsf{A}}}\equiv \left[ \begin{matrix} U & H & 0 & 0 & 0 \\ g(1-\beta) & U & \dfrac{g h}{\rho_w} & 0 & 0 \\ 0 & \varrho\chi & u & \varrho & 0 \\ g-\dfrac{\rm \pi}{\varrho h} & 0 & 0 & u & \dfrac{1}{\varrho}\\ 0 & \varrho\chi\mathscr{C}^2 & 0 & \varrho\mathscr{C}^2 & u \end{matrix} \right],\end{equation}

where $\beta \equiv \varrho /\rho _w$ and $\chi \equiv H/h$ are, respectively, the air–water density ratio and ocean-atmosphere depth ratio.

4.2. Linear waves

Let $\boldsymbol {\mu }_0 \equiv [(\eta _1)_0,U_0,\varrho _0,u_0,{\rm \pi} _0]^{{\rm T}}$ denote a stationary solution to the quasi-linear equation (4.1), and $\boldsymbol {\mu }_1 \equiv [(\eta _1)_1,U_1,\varrho _1,u_1,{\rm \pi} _1]^{{\rm T}}$ be a perturbation around the stationary solution, such that $\boldsymbol {\mu } = \boldsymbol {\mu }_0 + \boldsymbol {\mu }'$, where $\boldsymbol {\mu }'\equiv \zeta \boldsymbol {\mu }_1$ denotes the scaled perturbations and $\zeta$ is a smallness parameter ($\zeta \ll 1$). Replacing $\boldsymbol {\mu } = \boldsymbol {\mu }_0 + \zeta \boldsymbol {\mu }_1$ in (4.1) and collecting terms based on their $\zeta$ order results in

(4.2)\begin{align} \zeta^0:\qquad\qquad\qquad\qquad\qquad\ \qquad\qquad\qquad\qquad\qquad {\boldsymbol{\mathsf{A}}}_0\frac{\partial\boldsymbol{\mu}_0}{\partial s} &=\boldsymbol{0},\end{align}

(4.3)\begin{align} \zeta^1:\qquad\qquad\qquad\qquad\qquad\qquad\quad\frac{\partial\boldsymbol{\mu}_1}{\partial t}+{\boldsymbol{\mathsf{A}}}_1\frac{\partial\boldsymbol{\mu}_0}{\partial s}+{\boldsymbol{\mathsf{A}}}_0\frac{\partial\boldsymbol{\mu}_1}{\partial s} &=\boldsymbol{0}, \end{align}

where ${\boldsymbol{\mathsf{A}}}_0$ and ${\boldsymbol{\mathsf{A}}}_1$ indicate that ${\boldsymbol{\mathsf{A}}}$ is evaluated using $\boldsymbol {\mu }_0$ and $\boldsymbol {\mu }_1$, respectively. Assuming a uniform base flow, i.e. $\partial \boldsymbol {\mu }_0/\partial s = \boldsymbol {0}$, (4.2) is satisfied and (4.3) simplifies to

(4.4)\begin{equation} \frac{\partial\boldsymbol{\mu}_1}{\partial t}+{\boldsymbol{\mathsf{A}}}_0\frac{\partial\boldsymbol{\mu}_1}{\partial s}=\boldsymbol{0}. \end{equation}

Plane-wave solutions to (4.4) are sought in the form $\boldsymbol {\mu }_1=\hat {\boldsymbol {\mu }} \exp [\mathrm {i}(\omega t-ks)]$, where $k\in \mathbb {R}_{>0}$ is the wavenumber, $\omega \in \mathbb {R}_{>0}$ is the angular frequency, $\mathrm {i}^2=-1$, and $\hat {\boldsymbol {\mu }}=[\hat {\eta }_1,\hat {U},\hat {\varrho },\hat {u},\hat {{\rm \pi} }]^{{\rm T}}$ is a vector of eigenfunctions for the primitive variables. Substituting $\boldsymbol {\mu }_1$ using the eigenfunctions in (4.4) yields

(4.5)\begin{equation} ({\boldsymbol{\mathsf{A}}}_0 - \lambda{\boldsymbol{\mathsf{I}}}_5)\hat{\boldsymbol{\mu}} = \boldsymbol{0}, \end{equation}

with $\lambda \equiv \omega /k$. Non-trivial solutions to (4.5) exist if

(4.6)\begin{equation} \det{({\boldsymbol{\mathsf{A}}}_0 - \lambda{\boldsymbol{\mathsf{I}}}_5)}= \lambda (-\lambda^4 + (\mathscr{C}_w+\mathscr{C}_0^2)\lambda^2 + \mathscr{C}_w[\mathscr{C}_0^2({\beta_0}-1)+{\beta_0} \mathscr{C}_a - \mathscr{C}_T])=0, \end{equation}

where quiescent base flows for air ($U_0=0$) and water ($u_0=0$) have been assumed, $\mathscr {C}_w\equiv \sqrt {g H_0}$ and $\mathscr {C}_a\equiv \sqrt {g h_0}$ are the base flow gravity-wave speeds in the ocean and atmosphere layers, respectively, and $\mathscr {C}_0 \equiv \sqrt {\gamma {\rm \pi}_0/\varrho _0}$ and $\mathscr {C}_T\equiv \sqrt {{\rm \pi} _0/\varrho _0}$ are the base flow sound speed and Newton's sound speed in the atmospheric layer and $\beta _0 = \varrho _0/\rho _w$. The roots of (4.6) are

(4.7)\begin{equation} \lambda=0\quad \text{and} \quad \lambda^2 = \tfrac{1}{2} \left( \mathscr{C}_w^2 + \mathscr{C}_0^2 \pm \sqrt{ (\mathscr{C}_w^2 - \mathscr{C}_0^2)^2 + 4{\beta_0} \mathscr{C}_w^2 (\mathscr{C}_0^2 - \mathscr{C}_T^2 + \mathscr{C}_a^2) } \right), \end{equation}

giving rise to five eigenvectors. Each eigenvalue–eigenvector pair forms an eigenmode, namely,

(4.8) \begin{align} \left. \begin{array}{@{}l@{}} \mathscr{T} = \left\{ (\lambda,\hat{\boldsymbol{\mu}})\in\mathbb{R} \times \mathbb{R}^5:\lambda=0, \hat{\boldsymbol{\mu}} = \alpha \left[\begin{matrix} h_0 \\ 0 \\ \rho_w({\beta_0}-1) \\ 0 \\ {\rm \pi}_0-\varrho_0g h_0 \end{matrix}\right]\right\},\\ \mathscr{A^\pm} = \left\{ (\lambda,\hat{\boldsymbol{\mu}})\in\mathbb{R} \times \mathbb{R}^5: \lambda={\pm}\sqrt{\dfrac{1}{2}(X + \sqrt{Y^2+Z} )}, \hat{\boldsymbol{\mu}} = \alpha \left[\begin{matrix} \varphi H_0 \\ \varphi \lambda \\ \varrho_0 \left(\dfrac{\lambda^2}{\mathscr{C}_0^2}+\dfrac{\mathscr{C}_w^2}{\mathscr{C}_0^2}({\beta_0}-1)\right) \\ \lambda \left(\dfrac{\lambda^2}{\mathscr{C}_0^2}-\dfrac{\mathscr{C}_w^2}{\mathscr{C}_0^2}\right) \\ \varrho_0\mathscr{C}_0^2\left(\dfrac{\lambda^2}{\mathscr{C}_0^2}+\dfrac{\mathscr{C}_w^2}{\mathscr{C}_0^2}({\beta_0}-1)\right) \end{matrix}\right]\right\}, \\ \mathscr{G^\pm} = \left\{ (\lambda,\hat{\boldsymbol{\mu}})\in\mathbb{R} \times \mathbb{R}^5: \lambda={\pm}\sqrt{\dfrac{1}{2}(X - \sqrt{Y^2+Z} )},\ \hat{\boldsymbol{\mu}} = \alpha \left[\begin{matrix} \varphi H_0 \\ \varphi \lambda \\ \varrho_0 \left(\dfrac{\lambda^2}{\mathscr{C}_0^2}+\dfrac{\mathscr{C}_w^2}{\mathscr{C}_0^2}({\beta_0}-1)\right) \\ \lambda \left(\dfrac{\lambda^2}{\mathscr{C}_0^2}-\dfrac{\mathscr{C}_w^2}{\mathscr{C}_0^2}\right) \\ \varrho_0\mathscr{C}_0^2\left(\dfrac{\lambda^2}{\mathscr{C}_0^2}+\dfrac{\mathscr{C}_w^2}{\mathscr{C}_0^2}({\beta_0}-1)\right) \end{matrix}\right]\right\}, \end{array}\right\} \end{align}

where $X \equiv \mathscr {C}_w^2 + \mathscr {C}_0^2$, $Y \equiv \mathscr {C}_w^2 - \mathscr {C}_0^2$, $Z \equiv 4{\beta _0} \mathscr {C}_w^2(\mathscr {C}_0^2-\mathscr {C}_T^2+\mathscr {C}_a^2)$, $\alpha$ is any non-zero real number and $\varphi \equiv {\beta _0} \mathscr {C}_a^2/\mathscr {C}_0^2$.

Mode $\mathscr {T}$ is a stationary mode (relative to the air flow) reflecting the transport of thermal waves by the wind, characterised by the depth-averaged pressure ($\hat {{\rm \pi} }$) and density ($\hat {\varrho }$) fluctuations at constant Favre-averaged temperature ($\tilde {T}$). In the ocean, this mode induces a perturbation of the water column in phase opposition with the atmospheric pressure fluctuation in the standard atmosphere. Modes $\mathscr {A}^\pm$ are two acoustic-like modes propagating in opposite directions. In the limit of a zero-depth ocean, they converge to true acoustic modes in the atmosphere, evaluated using the layer-averaged density and pressure, i.e. $\lim _{H_0\to 0}(\lambda,\hat {\boldsymbol {\mu }}) = (\pm \mathscr {C}_0,\alpha [0,\pm \varphi \mathscr {C}_0,\varrho _0,\pm \mathscr {C}_0,\varrho _0\mathscr {C}_0^2]^{\textrm {T}})$ (and the gravity modes, discussed next, vanish). Modes $\mathscr {G}^\pm$ are two gravity-like modes propagating in opposite directions. In the limit of a zero-thickness atmosphere, they converge to the classical gravity waves in the ocean, i.e. $\lim _{h_0\to 0}(\lambda,\hat {\boldsymbol {\mu }}) = (\pm \mathscr {C}_w,\alpha [H_0,\pm \mathscr {C}_w,0,0,0]^{\textrm {T}})$ (and the acoustic modes vanish). These ‘gravito-acoustic’ modes (in the spirit of magneto-acoustic modes in plasmas), not to be confused with vertically dependent ‘acoustic-gravity waves’ discussed in, e.g. Vallis (Reference Vallis2017), are commented upon further in the following section in the context of planet Earth.

For comparison, the eigenmodes of the OWC model are derived here. The system in (1.1), augmented with the wave propagation equation for the ground/sea-level pressure signal, can be recast in a one-dimensional curvilinear system (with curvilinear coordinate $s$) and linearised around a base flow at rest with a uniform sea floor yielding

(4.9)\begin{gather} \frac{\partial\boldsymbol{q}_1}{\partial t} + {\boldsymbol{\mathsf{A}}}_{owc} \frac{\partial\boldsymbol{q}_1}{\partial s} = \boldsymbol{0}, \quad\mbox{with:}\ \boldsymbol{q}_1\equiv\left[ \begin{matrix} {(}\eta_1{)_{1}}\\ U_1\\ v_1\\ p_1 \end{matrix} \right], \quad{\boldsymbol{\mathsf{A}}}_{owc} \equiv\left[ \begin{matrix} 0 & H_0 & 0 & 0 \\ g & 0 & 0 & 1/\rho_w\\ 0 & 0 & 0 & 1/\rho_a\\ 0 & 0 & \rho_a c_a^2 & 0 \end{matrix}\right],\end{gather}

where, respectively, $U_1$ and $v_1$ are the velocity perturbations in the ocean and at ground level in the atmosphere, $p_1$ and $c_a$ are the ground level atmospheric pressure perturbation and its fixed propagation speed, $\rho _a$ and $\rho _w$ are the ground level air and water density, and ${(}\eta _1{)_{1}}$ is the sea-surface perturbation around the constant water depth $H_0$. Seeking plane-wave solutions $\boldsymbol {q}_1 = \hat {\boldsymbol {q}}\exp [\mathrm {i}(\omega t - k s)]$ to (4.9), where $(\omega,k) \in \mathbb {R}^2$ are, respectively, the angular frequency and wavenumber, reveals that the OWC system propagates two pairs of eigenmodes, defined as

(4.10)\begin{gather} \mathscr{A}^\pm_{owc} = \left\{(\lambda,\hat{\boldsymbol{q}})\in\mathbb{R} \times \mathbb{R}^4:\ \lambda={\pm} c_a,\ \hat{\boldsymbol{q}} =\alpha\left[\begin{matrix}\xi H_0\\\pm\xi c_a \\ \pm c_a \\ \rho_a c_a^2\end{matrix}\right]\right\}, \end{gather}

(4.11)\begin{gather}\mathscr{G}^\pm_{owc} = \left\{(\lambda,\hat{\boldsymbol{q}})\in\mathbb{R} \times \mathbb{R}^4:\ \lambda={\pm} \mathscr{C}_w,\ \hat{\boldsymbol{q}}=\beta\left[\begin{matrix}H_0\\ \pm \mathscr{C}_w\\ 0\\ 0 \end{matrix}\right]\right\}, \end{gather}

where $\lambda \equiv \omega /k$, ($\alpha$, $\beta$) $\in \mathbb {R}_{\neq 0}^2$, $\mathscr {C}_w \equiv \sqrt {g H_0}$ (retaining its definition as in the TWC case) and $\xi \equiv (\rho _a/\rho _w)c_a^2/(c_a^2-\mathscr {C}_w^2)$. The pair $\mathscr {A}^\pm _{owc}$ corresponds to left- ($-$) and right- ($+$) going fluctuations in both the atmosphere and ocean at the set speed $c_a$, whilst the pair $\mathscr {G}^\pm _{owc}$ corresponds to the classical (tsunami) left-/right-going gravity modes in the ocean, with no associated atmospheric fluctuations.

5.1. Eigenmodes in a standard atmosphere

The standard density, pressure and temperature profiles of Earth's atmosphere are shown in figure 4, together with the depth-averaged density ${\varrho }_0^*$, pressure ${\rm \pi} _0^*$ and the Favre-averaged temperature ($\tilde {T}_0^*$) for $h^*_0$ ranging from ground level to $80\ \mbox {km}$. Here, the superscript $^*$ denotes dimensional quantities. Materials from the Tonga eruption are reported to have reached heights in excess of $50$ km (Carr et al. Reference Carr, Horváth, Wu and Friberg2022) and signatures from higher layers (e.g. GNSS measurements of ionospheric disturbances Wright et al. Reference Wright2022) suggested that a layer at least $80$ km thick was perturbed. Remarkably, the Favre-averaged temperature is found to remain nearly constant for heights beyond $20\ \mbox {km}$ (see figure 4), which means that $\mathscr {C}_0^*$ becomes nearly constant for $h^*_0 > 20\ \mbox {km}$. This value is ${\sim } 317\ \mbox {m}\ {\cdot }\ \mbox {s}^{-1}$, in good agreement with the observed Lamb-wave propagation speed of $318 \pm 6\ \mbox{m}\ {\cdot }\ \mbox {s}^{-1}$ (Wright et al. Reference Wright2022). Thus, considering the atmosphere to be $h_0^* \sim 80\ \mbox {km}$ thick seems appropriate for sub-ionospheric considerations for the Tonga event. Note that this choice is compatible with the choice in (3.10) as $p^*|_{\eta _2} \ll \varrho _0^* {g_o} h_0^*$. Considering $h_0^* = 80$ km and $\gamma = 1.4$, the logarithmic average values of the standard atmosphere are ${\varrho }_0^* = 0.129\ \mbox {kg}\ {\cdot }\ \mbox {m}^{-3}$, ${\rm \pi} _0^* = 9.30\times 10^3\ \mbox {Pa}$ and $\mathscr {C}_0^* = 317\ \mbox {m}{\cdot }\mbox {s}^{-1}$. Together with $\rho _o = 10^3\ \mbox {kg} {\cdot } \mbox {m}^{-3}$, $g_o = 9.81\ \mbox {m}{\cdot }\mbox {s}^{-2}$ and $\ell _o = 3668\ \mbox {m}$ (the average water depth on Earth), the eigenvalues for modes $\mathscr {A}^\pm$ and $\mathscr {G}^\pm$ are evaluated in figure 5.

Figure 4. (a) Vertical profiles of pressure ($p^*$: –, ${\rm \pi} ^*_0:$ -- ) and density ($\rho ^*$: - (red), $\varrho ^*_0$: -- (red)). (b) Vertical profiles of temperature (–), Favre-averaged temperature (--), local speed of sound (- (red)) and $\mathscr {C}_0^*$ (-- (red)). All values are computed using the international standard atmosphere.

Figure 5. Propagation speeds of modes $\mathscr {A}^+$ and $\mathscr {G}^+$ for the OWC (a) and TWC (b) models applied to the standard atmosphere as a function of the local water depth. The probability of observing the water depth on Earth is shown by the filled blue curve (labelled ‘pdf’, normalised by its highest value), where ‘CD’ is the Challenger Deep point (Mariana Trench). The wave speed for mode $\mathscr {A}^+$ in the OWC is arbitrarily set. Results for $280\ \mbox {m}\,\mbox {s}^{-1}$, $300\ \mbox {m}\,\mbox {s}^{-1}$ and $320\ \mbox {m}\,\mbox {s}^{-1}$ are shown, following the work of Kubota et al. (Reference Kubota, Saito and Nishida2022) (KSN). The bottom figures give some of the normalised eigenfunctions for the OWC (a,c) and TWC (b,d) models: (i) is any of $\hat {\varrho }/\varrho _0$, $\hat {u}/\mathscr {C}_0 \hat {{\rm \pi} }/(\varrho _0\mathscr {C}_0^2)$, $\hat {T}/((\gamma -1)T_0)$, (ii) is $\hat {\eta }/(\varphi H_m)$ for both the $\mathscr {A}^+$ and $\mathscr {G}^+$ modes, (iii) is $\hat {\eta }/H_m$ for the OWC $\mathscr {G}^+$ mode, where $H_m = 20\ \mbox {km}$ (the range of water depths shown here). Finally, the water-surface displacement $|\eta ^\prime |$ (in blue) associated with a ground pressure fluctuation of $2.5\ \mbox {hPa}$ from the $\mathscr {A}^+$ mode is shown for the OWC, TWC and hydrostatic (STA) models (see equations in § 5.2). The thin blue line labelled $\mbox {(*)}$ in the bottom right panel is the OWC result with a wave speed set to $\mathscr {C}_0$. See comments in the text.

Two distinct regimes emerge. Below a critical water depth, $H_c^*$, the wave speed of the acoustic mode is nearly independent of the water depth, whereas the wave speed of the gravity mode increases with the water depth (proportionally to $\sqrt {g_o H^*}$). Conversely, above the critical depth, the acoustic mode speeds up (proportionally to $\sqrt {g_o H^*}$) whereas the gravity mode saturates at a speed that never exceeds the lowest speed of the acoustic mode. This implies that at no point the gravity mode can be in resonance with the acoustic mode (this differs from the resonance that can be achieved with the considerations in Proudman (Reference Proudman1929), and subsequent OWC models, where the pressure-wave speed can match the gravity-wave speed exactly). The absence of a resonance is guaranteed by the fact that $Z$ is always positive.

The critical depth, $H_c$, corresponds to the $H_0$ value for which the wave-speed of $\mathscr {A}^\pm$ comes the closest to that of $\mathscr {G}^\pm$. By inspection of the acoustic- and gravity-mode eigenvalues, they approach each other if the inner square root is minimised, that is, when $f(H_0) = (g H_0 - \mathscr {C}_0^2)^2 + 4{\beta _0} g H_0 (\mathscr {C}_0^2 - \mathscr {C}_T^2 + \mathscr {C}_a^2)$ reaches a minimum. This occurs if $H_0 = \mathscr {C}_0^2(1-2{\beta _0})/g +2{\beta _0}(\mathscr {C}_T^2-\mathscr {C}_a^2)/g \equiv H_c$, which in dimensional form gives

(5.1)\begin{equation} H_c^* = \frac{{\mathscr{C}_0^*}^2}{g_o}\left(1-2\frac{\bar{\varrho}_0^*}{\rho_w^*}\right) + 2 \frac{\bar{\varrho}_0^*}{\rho_w^*g_o}({\mathscr{C}_a^*}^2-{\mathscr{C}_T^*}^2) \approx 10\ \mbox{km}. \end{equation}

Remarkably, this critical-depth value coincides with the deepest water bodies on Earth. This means that the vast majority of oceans on Earth will only experience the subcritical regime, where the acoustic mode propagates at a nearly constant speed $\mathscr {C}_0^* \approx 317\ \mbox {m}{\cdot }\mbox {s}^{-1}$ (as observed, e.g. Wright et al. Reference Wright2022).

Moving to the eigenvector components, the subcritical regime is dominated by the acoustic mode in the atmosphere, whereas the ocean's surface exhibits the footprint of the acoustic mode on the water-height disturbance in addition to the two (classical) gravity waves. If oceans were deep enough to be in the supercritical regime, the atmosphere would instead be dominated by fluctuations from the gravity modes (from the ocean) at a near-constant speed $\sqrt {g_o H_c^*} \approx 313\ \mbox {m}{\cdot }\mbox {s}^{-1}$, and oceans would experience surface waves from these gravity waves in addition to faster acoustic modes with almost no associated fluctuations in the atmosphere. The only places where this could be observed are above the deepest oceanic trenches.

5.2. $\mathscr {A}$-mode associated water-height fluctuations

The amplitude of the water-wave height, $|\eta ^\prime | =|\zeta \hat {\eta }|$, relative to the depth-averaged pressure-fluctuation amplitude, $|{\rm \pi} ^\prime |=|\zeta \hat {{\rm \pi} }|$, is of particular interest since it quantifies the relative importance of the potential energy in the ocean with that of the internal energy in the atmosphere. However, given that ground pressure fluctuations $p^\prime (\eta _1)$ are more readily measurable (and it is the term identified as performing the work on the sea surface in OWC models), for ease of interpretation, it is useful to relate it to depth-averaged pressure fluctuations. Linearising the pressure boundary condition (3.10) and using the relation between the $\mathscr {A}$-mode-related pressure and density fluctuations yields

(5.2)\begin{equation} |p'(\eta_1)| = \zeta | \hat{\varrho}g h_0 - \varrho_0g \hat{\eta} | = \frac{|{\rm \pi}^\prime|}{|\hat{\rm \pi}|} | \hat{\varrho}g h_0 - \varrho_0 g\hat{\eta} |. \end{equation}

Although the ratio of $|p^\prime (\eta _1)|$ over $|{\rm \pi} ^\prime |$ depends on the local water depth $H_0$ (since the eigenfunctions depend on $H_0$), the dependency is found to be weak when evaluated over Earth's bathymetry. The ratio for the standard atmosphere can be estimated as $|p^\prime (\eta _1)|/|{\rm \pi} ^\prime | \sim 7.8$ for any depth up to the critical depth.

Let $\delta p^*$ represent the ground/sea-level pressure fluctuation following the eruption and $\delta H^*$ be the associated sea-surface displacement. If using a purely hydrostatic argument, $\delta p^* \sim \rho _o g_o\,\delta H^*$, leading to the widespread result $\delta H^*/\delta p^* \sim 1\ \mbox {cm}{\cdot }\mbox {hPa}^{-1}$ (Röbke & Vött Reference Röbke and Vött2017), whereas the eigenmode analysis gives

(5.3)\begin{equation} \frac{\delta H^*}{\delta p^*} = \frac{\zeta |\hat{\eta}|}{|p^\prime(\eta_1)|}\frac{\ell_o}{\rho_o g_o\ell_o} = \frac{1}{\rho_o g_o}\frac{|\hat{\eta}|}{| \hat{\varrho}g h_0 - \varrho_0g \hat{\eta}|}.\end{equation}

This ratio is used to plot the resulting water-height fluctuation for a sea-level pressure perturbation of $2.5$ hPa over the range of depths on Earth in figure 5 (bottom right panel). Remarkably, the resonance observed in the OWC does not appear in the TWC as the $\delta H^*/\delta p^*$ ratio transitions continuously from the subcritical to the supercritical regime. At the critical depth, the eigenmodes are almost collinear and, therefore, assume similar properties, i.e. the distinction between acoustic and gravitational origin is more difficult. This translates to perturbations mostly on the water surface. This is interpreted as a directional energy transfer towards potential energy in the ocean. For Earth, $\delta H^*/\delta p^* \sim 3\times 10^1\ \mbox {cm}{\cdot }\mbox {hPa}^{-1}$ is observed at the critical height, i.e. one order of magnitude greater than the usual hydrostatic argument. In shallow oceans, $\delta H^*/\delta p^*$ decreases substantially and the atmospheric pressure wave is not associated with large water-height perturbation. This leads to very local air to water energy transfers. Whilst a ${\sim } 2.5 \ \mbox {hPa}$ pressure wave from the eruption hardly disturbs the ocean over continental shelves, it is found to travel with a $70\ \mbox {cm}$ high wave over the deepest oceanic trenches. Such wave heights are comparable with the more common tsunami generated by a sudden seafloor motion (Ward Reference Ward2002) and explain the significant impact of the eruption compared with what is predicted by hydrostatic considerations.

5.3. One-dimensional atmospheric wave propagation

Given the discussion above for an idealised atmosphere, its applicability to atmospheric wave propagation around the globe following the Tonga event is investigated. As illustrated in figure 1, the atmospheric wave initially presents a remarkably circular propagation away from the source. To better visualise this, the IR data of figure 1 is mapped to a cylindrical map projection whose poles lie at the Hunga Tonga-Hunga Ha'apai (HTHH) volcano and its antipode; thus, each point is identified by its forward azimuth and its distance from HTHH. The signal is then averaged at an iso-distance over azimuth segments to increase the signal-to-noise ratio and allow the wave position to be more clearly identified and manually extracted (see Appendix C for details). The result is shown at various times on 15 January 2022 in figure 6 (in this projection, a rectilinear wave corresponds to a circular wave on a globe). Within the first approximately 4 hrs after the explosion, the wave propagation is highly symmetric (a maximum difference of ${\sim } 2 \times 10^2\ \mbox {km}$ is observed when comparing the extracted wave position to a perfectly circular wave at 8:00UTC). However, at later times, the wavefront exhibits significant departure from a circular propagation as its starts to kink (by 19:00UTC, the maximum deviation from circular propagation reaches ${\sim } 10^3\ \mbox {km}$). Two first-order effects a priori responsible for the wave deformation are identified here: wave propagation speed modification due to underlying thermodynamic/topographic variations and that due to underlying atmospheric flows.

Figure 6. Projection of the sector-averaged global IR data (9.6 $\mathrm {\mu }$m band) onto a cylindrical map with poles at HTHH and its antipode, respectively, at times 7:00, 11:30 and 15:30 (panels (b)–(d), ordered from top to bottom). The extracted wave position is shown in cyan. For clarity, the IR data are faded away from the extracted wave position. Topography is shown in black with the main mountain ranges in colour. Panel (b) contains region indications (AU: Australia, SA: South America, EU: Europe, AF: Africa, AS: Western Asia). A 10$^\circ$ silver patch is added around the North (denoted by ‘N’) and South poles. A quiver plot of velocity projected on the isoazimuth lines (see Appendix D for magnitudes) is shown in white in the panel (a).

To evaluate the impact of each effect, realistic data from the day (ERA5 data that is a fusion of observation and simulation data, see, e.g. Hersbach et al. Reference Hersbach2018) is used to integrate the equation governing wave position in time along one-dimensional lines from HTHH to its antipode (see Appendix D for details). The predicted wave positions are shown in figure 7 where, for reference, the result using uniform atmospheric conditions (using international standard atmosphere values) with a quiescent flow is also shown. Taking into account solely thermodynamic and topographic variations highlights faster wave propagation closer to the equator (where the Favre-averaged temperature is higher) and a slower propagation near the poles (where the Favre-averaged temperature is lower) as well as more subtle effects due to the highest mountain ranges (e.g. Himalayas and Andes). However, these effects alone are not enough to account for much of the wave deformation. Taking into account the underlying atmospheric flows highlights the impact of the main features of atmospheric planetary circulation: circulation around each of the poles is predominantly co (counter) current with the atmospheric wave propagation as along positive (negative) azimuth lines resulting in locally advanced portions of the waves (centred around azimuths ${\sim } 50^{\circ }$ and ${\sim } 150^{\circ }$) as well as locally retarded portions (centred around azimuths ${\sim } -40^{\circ }$ and ${\sim } -160^{\circ }$). The result of the combined effect of both contributions is illustrated in the last panel of figure 7. The start time is set at 4:30 UTC that was obtained by computing the time shift necessary to minimise the $L_2$ norm between the extracted wave position and the integrated solution at 10:00 UTC. This value is in agreement with the estimated initial time of the main explosion (Wright et al. Reference Wright2022).

Figure 7. Integration of (D1) with a uniform atmosphere and stationary flow (a), a non-uniform atmosphere and stationary flow (b), a uniform atmosphere and non-uniform velocity field (c) and non-uniform atmosphere and velocity field (d) shown in grey lines. Atmospheric values and velocity fields are set using ERA5 values from 15 January 2022. The extracted wave position is shown in black (with portions of medium and low confidence in the extraction shown in orange and red dashed lines, respectively). Note that as time goes on, due to the loss of circular coherence, the signal-to-noise ratio decreases and the wave is harder to extract, thus, larger portions of the curve have a lower confidence rating. The comparison is shown at every hour from 05:00UTC to 19:00UTC.

Remaining deviations (e.g. the excessive wave lag around azimuth ${\sim } -140^{\circ }$) are attributed to multiple factors not taken into account in the current integration, e.g. humidity and atmospheric composition (modifying the local speed of sound), lack of precision of the global data for the 15 January 2022 and multi-dimensional effects. For the latter, see the results of 2-D simulations in § 5.4.

5.4. Two-dimensional numerical simulations of the event

With the aim of illustrating how 2-D effects, spherical geometry and realistic bathymetry affect the one-dimensional considerations presented above, global 2-D simulations taking into account the non-uniform atmosphere and sea floor are presented in the following subsections.

5.4.2. Base flow and source condition

Given the importance of the underlying thermodynamic quantities and velocity fields observed when integrating the eigenmode propagation, the same 2-D base flows are taken into account in the simulation. In addition, realistic high-resolution bathymetry and topography from the general bathymetric chart of the oceans (GEBCO) are used to represent the ocean floor and terrain height. With the intention of preventing unrealistic wave generation and propagation over the North Pole, the bathymetry data are supplemented with ice coverage for both poles for January 2022 that was obtained from the National Snow and Ice Data Center. The governing equations advanced in time are of the form

(5.4)\begin{equation} \frac{\partial \boldsymbol{\mu}}{\partial t} = {\text{RHS}}(\boldsymbol{\mu}). \end{equation}

However, the underlying initial fields of $\boldsymbol {\mu } (t=t_i^- )$, where $t_i^-$ is the initial time prior to adding the explosion source, are not a steady-state solution to the governing equations. Thus, the right-hand side (RHS) resulting from these quantities, ${\textbf {RHS}}(\boldsymbol {\mu }(t=t_i^-))$, prior to the addition of the perturbation modelling the explosion, is computed at the start of the simulation and then subtracted from the right-hand side for the rest of the time integration to ‘freeze’ the base flow. Such a strategy is employed in, e.g. linear stability analysis to study the response of systems where an imposed time-averaged field (which is not a steady-state solution to the governing equations) is disturbed by small-amplitude fluctuations (see, e.g. Touber & Sandham Reference Touber and Sandham2009).

The aim of these simulations is not to reproduce the short-time nonlinear shock-wave dynamics associated with the initial explosion. Attention is focused on the time when the atmospheric wave has reached a quasi-linear behaviour and established its vertical profile. To reflect this, the volcano source is given a spatial and temporal support: the spatial scale $\varLambda$ is chosen to provide one scale separation from when the wave is still a shock wave and the temporal scale $\tau$ is obtained from numerous ground-based pressure measurements. The shape of the spatial support is given by a 2-D Gaussian of standard deviation $\varLambda$. This is compatible with the initially symmetric nature of the atmospheric wave observed in figure 1. The shape of the temporal support is guided by the observed ground pressure signal (e.g. those of figure 12) and the following criteria are retained: the slope and value of the function at $t=0$ and $t=\tau$ must be zero and the function must reach an amplitude $p^+$ at $t=\tau /4$ and $p^-$ at $t=3\tau /4$. A fifth-order polynomial is constructed to satisfy these constraints. Shock propagation distance estimates (e.g. Lynett et al. Reference Lynett2022) lead to the choice of $\varLambda = 50\ \mbox {km}$ and ground pressure measurement processing (e.g. those used in figure 2) lead to the choice of $\tau =34\ \mbox {min}$, $p^+= 5.2\ \mbox {hPa}$ and $p^- = -0.1 p^+$.

5.4.4. Results

5.4.4.1 Atmospheric wave propagation

As illustrated in figure 1, despite its remarkably circular initial shape, IR satellite data clearly shows the wave deforming as it travels to Tonga's antipode. Integration of the eigenmodes with realistic conditions on the day, i.e. with underlying thermodynamic and velocity fields (see figure 7), motivated the investigation of 2-D effects. Figure 8 presents a comparison between the brightness temperature fluctuations (observed by geostationary satellites) and the ${\rm \pi} '$ field obtained from the numerical simulation at 1 hr intervals from 15:30 to 18:30 UTC. As the wave propagates around the world, its interactions with the topography and the base flow cause the initially circular front to distort. These distortions are remarkably well captured by the simulation. For example, as observed in one dimension, the passage near the South Pole locally advances/retards portions of the wave such that, in two dimensions, a wave superposition is observed by the time the wave reaches the west coast of Namibia and Angola leading to a local maximum of the fluctuation amplitudes. A similar behaviour is observed for the portion of the wave off the northwest coast of Africa. This can explain why the wave is more clearly visible at the corresponding locations in the IR data as there will be an associated local extremum of $\tilde {T}'$. In addition, the interaction of the wave with the topography is apparent in the simulation data where reflected waves from the passage of the main wave on the Andes mountain range can be seen propagating to the west (second pressure level fluctuations in figure 8). This was also captured by a geostationary satellite and illustrated in, e.g. Wright et al. (Reference Wright2022) who processed GOES-16 data to make the reflected wave visible. This reflected wave goes on to interact with the sea surface that was already disturbed by the passage of the primary wave. As evidenced by the plot, many other smaller amplitude reflected waves were generated by interactions with other topographic features (e.g. peaks in Hawaii, New Zealand) and propagate throughout the Pacific. Supplementary movies 1–4 (available at https://doi.org/10.1017/jfm.2023.131) provide dynamical views of the simulation results.

Figure 8. Brightness temperature (from the $9.6\ \mathrm {\mu }\mbox {m}$ centred spectral band) fluctuations ($\mathcal {T}^\prime$) (a,c,e,g) from geostationary satellites (see Appendix A for details), low-pass filtered with a $75\ \mbox {km}$ cutoff and saturated at $\pm 0.01\ \mbox {K}$ for the purpose of highlighting the thermal-wave position on the world map. Colours give the topography with the same scale as in figure 1. Time is given in UTC on 15 January 2022. Depth-averaged pressure-fluctuation fields (${\rm \pi} '$) from the numerical simulation of the TWC model at the same instants as the thermal wave are provided (b,d,f,h). Two levels are shown: $\pm 0.2$ hPa in light/dark and $\pm 0.01$ hPa in blue/yellow. For the latter, the field is clipped to 10 500 km from Tonga so as to not obscure the first field. The location of selected ground pressure measurements in figure 12 are illustrated by the coloured dots. The colour indicates which group of sensors the measurement belongs to (green: Sensor Com. Archive, red: Weather News Inc, blue: NOAA/ASOS). The snapshots correspond to times when the main wave is focusing towards HTHH's antipode in the Sahara (yellow marker). The TWC model correctly predicts the wave location, including the pinches along its front. The trailing waves in the satellite data correspond to subsequent eruptions following the most energetic one, see, e.g. Wright et al. (Reference Wright2022).

5.4.4.2 Historical maximum of absolute water displacement

During the numerical simulations, the historical maximum of the atmosphere–water surface interface perturbation, $|\eta _1'|$, is updated at regular intervals in time (every 30 physical seconds) over the whole domain. Figures 9 and 10 propose a comparison of this historical maximum 18 hrs after the explosion between three different scenarios: a case where only the water layer is advanced in time (ZWC model) with a source contribution due only to the $\eta _1$ contribution of the $\mathscr {A}^{\pm }$ mode over the temporal support; a case where the OWC model is integrated in time according to the source conditions of § 5.4.3, and a case where the TWC model is integrated in time according to the source conditions of 5.4.2. Comparing the ZWC scenario to the other two highlights the role of the interaction of the atmospheric wave with the non-uniform water depth. Indeed, the energy injected into the water layer at the source during the perturbation time scale $\tau$ by the $\mathscr {A}^{\pm }$ modes cannot account alone for the observed wave heights or propagation times (note that this is not equivalent to a comment on any $\mathscr {G}$-mode contributions, see discussion below regarding energy injected per layer considerations). When the atmospheric forcing is present, significant $|\eta _1|$ fluctuations (i.e. > 4 cm) can be broken down into three categories.

1. (i) As an $\mathscr {A}$ mode travels over a deep basin with a uniform depth, it propagates with it an $\eta _1$ fluctuation proportional to its associated pressure-fluctuation amplitude according to the relation graphed in figure 5 (this signal is observed in sensor data, discussed next).

2. (ii) The $\mathscr {A}$ mode transfers part of its energy to $\mathscr {G}$ modes as it interacts with steep changes in bathymetry. This is particularly clear near the volcano as the atmospheric wave travels over the Tonga and Kermadec trenches: they are the loci of generation of large $\mathscr {G}$-mode associated $\eta _1$ fluctuations that, once seeded, then propagate to the southwest towards South America's eastern coast and Antarctica. The refraction problem, i.e. the theory governing the mode-to-mode energy transfer at step changes in water depth, is derived below.

3. (iii) The local constructive superposition of multiple $\mathscr {G}$-mode associated $\eta _1$ fluctuations.

Figure 9. (a) $\mathscr {A}$ associated $\eta _1'$ for a 2.5 hPa sea-level pressure perturbation at every submerged point on the globe (various deep basins are numbered with ‘BX’ notation and referred to in the text). (b) Historical maximum of the absolute water-height fluctuation over the 18 hrs following the explosion for the TWC model. The cyan dots represent the locations of the DART sensors used in figure 13.

Figure 10. Historical maximum of the absolute water-height fluctuation over the 18 hrs following the explosion for (a) the ZWC model and (b) the OWC model. Note the different scales between the two plots. The cyan dots represent the locations of the DART sensors used in figure 13.

Figure 9 proposes a comparison between the $\mathscr {A}$-mode associated $|\eta _1'|$ for a uniform global ground pressure perturbation of 2.5 hPa and simulation results over the same region 18 hrs after the explosion to better discern which of the aforementioned effects is responsible for local maxima in the historical maximum maps. Note that due to the spherical shell geometry and deformation of the wave over time, the local pressure perturbation evolves in space, thus, an equal depth region close and far from Tonga will not see the same $|\eta _1'|$ (unlike what is assumed in the top panel of figure 9). The map at 18 hrs is the superposition of all the effects above; figure 11 and supplementary movies 1–4 are provided to show how the map is formed in time over selected areas around the globe. Nevertheless, regions where one effect dominates over the others can be identified (numbered B notation refers to those used in figure 9) as follows.

1. (a) Prominent examples of (i) can be seen after the atmospheric wave crosses the South Pole and begins to pinch, locally inducing a higher ${\rm \pi} '$; this leads to the streak of local maxima between Antarctica and South Africa over B7 (see figure 11 row 4). A similar phenomenon can be observed as the atmospheric wave travels between the west coast of Brasil and Senegal. In addition, the deep basin B4 (northwest pacific basin), off the west coast of Japan, registers large $\mathscr {A}$-mode associated $|\eta _1'|$ (see figure 11 row 2).

2. (b) The most prominent example of (ii) is at the Tonga–Kermadec trench (descent into B5), however, other notable examples can be found, e.g. at the step down from the coast of Argentina (into B6, see figure 11 row 4), the step down off the coast of Florida to the Atlantic Ocean (into B9, see figure 11 row 4), step down from the coast of Southern Australia towards Antarctica (into B2, see figure 11 row 2).

3. (c) Prominent examples of (iii) can be observed, e.g. in the region between Southern Australia and Antarctica as two $\mathscr {G}$ modes are generated at an angle by the passage of the $\mathscr {A}$ mode and interact to form a local maximum as they travel to the west (see figure 11 row 2). Other such interactions create the ray-like features in the historical maximum water-height disturbance map as $\mathscr {G}$ modes interact with Hawaii and the islands between it and the west coast of North America (see figure 11 row 3).

Figure 11. Time lapses of the simulation results for the TWC model, extracted from supplementary movies 1–4. Blue–red colours (see colourbar) show the maximum absolute sea-surface displacement measured from the initial (at rest) condition up to the time displayed on the top-left corner (given in UTC for 15 January 2022). The yellow/orange bands show the depth-averaged pressure fluctuation in the atmosphere where it exceeds $+5\ \mbox {Pa}$ (yellow) or is below $-5\ \mbox {Pa}$ (orange). The greyscale highlights the background topography, where both high mountain ranges and deep basins appear brighter. Shading effect is also applied to the topography and the instantaneous ocean waves to give depth effects. Each row comes from different camera paths along a constant azimuth from Tonga. See the movies for the full path from Tonga to its antipode. These views are used to illustrate the physical mechanisms at play discussed in the main text.

The effects in (i) and (ii) are idealised as they both consider a water depth that is piecewise uniform, but they can be used to comment on the more pronounced features of the historical maximum water-height disturbance maps. In practice, any variation of the water depth will be responsible for some amount of energy transferred between the layers. This is observed in the simulation data (figure 11 and supplementary movies 1–4) as the atmospheric wave generates many local sub-centimetre $|\eta _1|$ fluctuations as it propagates over the water away from large gradients of depth. Constructive interference between the various propagating waves can result in local maxima that will be highly sensitive to the source condition (e.g. when attempting to include a ‘more realistic’ frequency content in the source condition).

When comparing the OWC and TWC coupled results, qualitative similarities exist in the region near the initial perturbation where large $\eta _1$ fluctuations are generated at the Tonga–Kermadec trench and various ray structures are sent toward the northeast. Quantitatively, the models diverge as the OWC model shows greater $|\eta _1'|$ over larger regions than the TWC model. This is due, in part, to the form of the source condition that imposes a pure pressure perturbation to model the explosion that is not equivalent to imposing an $\mathscr {A}_{owc}$-type perturbation. Indeed, decomposing a pure pressure fluctuation into its eigenmode contributions shows that it is equivalent to imposing both $\mathscr {A}^\pm _{owc}$ and $\mathscr {G}^\pm _{owc}$ modes,

(5.5)\begin{equation} \begin{bmatrix} 0, 0, 0, 2 \alpha_f \rho_a c_a^2 \end{bmatrix}^{{\rm T}} = \left( \hat{\boldsymbol{q}}_{\mathscr{A}^+_{owc}} + \hat{\boldsymbol{q}}_{\mathscr{A}^-_{owc}} - \dfrac{\xi \alpha_f}{\beta} ( \hat{\boldsymbol{q}}_{\mathscr{G}^+_{owc}} + \hat{\boldsymbol{q}}_{\mathscr{G}^-_{owc}} ) \right), \end{equation}

where the eigenvector notations of § 4.2 have been used and $\alpha _f$ is a constant set based on the desired amplitude of the pressure perturbation.

The models further diverge as the atmospheric wave propagates away from the Pacific: unlike the uniform-amplitude circular wave imposed by the OWC model, the variations in the shape and amplitude of the atmospheric wave (in the TWC model) modify the location and intensity of large $|\eta _1|$ fluctuations. This is particularly noticeable in figure 9 in the southern Atlantic Ocean. As shown in figure 8, the amplitude of the atmospheric wave is weaker off the east coast of Brazil than it is closer to Africa, which leads to the overestimated (underestimated) $|\eta _1|$ fluctuations close (further) to Brazil and underestimated fluctuations along the southern and east coasts of Africa. Note also that a modification of the atmospheric wave shape will modify the angle of incidence when approaching step changes in water depth and, thus, modify the refraction properties. These aspects highlight the importance of accounting for atmospheric non-uniformities when attempting to provide approximately 24 hr time-scale predictions.

5.4.4.3 Sensor comparison

Two sets of ‘ground truth’ sensor comparisons are presented here. First, selected signals from an array of ground pressure sensors are shown in figure 12. For comparison, the ground pressure is extracted from dNami by computing the fluctuations of $p (\eta _1)$ according to (3.10). The signals illustrate how the proposed model is able to capture some of the subtle changes in the amplitude and shape of the atmospheric wave, e.g. as the wave approaches and passes over the United States, it is affected by a rapidly spatially varying $\boldsymbol {U}$ on the west United States coast and the Rocky mountain range resulting in a lower amplitude along the central latitude of the contiguous United States. Note how the simulated propagation speed of the wave, a consequence of the prescribed thermodynamic base flow and velocity fields, is in line with the observed propagation speed from the ground pressure signals. Secondly, signals from ocean bottom pressure sensors from the deep-ocean assessment and reporting of tsunamis (DART) network are gathered and presented in figure 13. These sensors are located in deep water (between 1.8 to 5.9 km) where a quantitative comparison between the data and the long-wave theory can be made (i.e. away from any coast-related wave steepening effects). For most of the signals, three regimes/arrival times have been identified (Lynett et al. Reference Lynett2022; Matoza Reference Matoza2022): (i) the arrival of the $\mathscr {A}$ mode, (ii) the arrival of the ‘principal’ $\mathscr {G}$ mode and (iii) secondary fluctuations generated by $\mathscr {A}$- and $\mathscr {G}$-mode interactions with the bathymetry. Note that the passage of the $\mathscr {A}$ mode does not always trigger the high-resolution recording mode of the DART sensors, therefore, some signals are less temporally resolved than others. For comparison, the ocean bottom pressure is extracted from dNami, for both OWC and TWC models, by computing the fluctuations of $p (\eta _0)$ according to $p(\eta _0) = p(\eta _1) + \rho _w g (\eta _1 - \eta _0)$. In both cases, the $\mathscr {A}$-mode-related pressure fluctuations are well predicted. As expected from the eigenmode analysis of figure 5, when running the OWC with the appropriate wave propagation speed, energy injection from the atmospheric layer into the water layer is similar to that of the TWC model. Note however that after the passage of the $\mathscr {A}$ mode, the OWC yields slightly greater $p(\eta _0)$ fluctuations than the TWC. This is due, in part, to the choice of source forcing in the OWC case discussed above. The inclusion of a correctly modelled source $\mathscr {G}$-mode contribution is discussed next.

Figure 12. Select ground pressure signals from sensors around the world (black) and signal extracted from dNami (red) within the first 18 hrs since the explosion. For each signal, a coloured dot indicates which group of sensors the measurement belongs to (green: Sensor Com. Archive, red: Weather News Inc, blue: NOAA/ASOS) with their locations displayed in figure 8 using the same colour coding. The azimuth from the sensor to the volcano is given in the bottom right-hand corner of each panel. The panels are sorted by increasing distance from the volcano from top left to bottom right.

Figure 13. Ocean bottom pressure from DART (black), dNami TWC (red), dNami OWC (blue). The DART code and sensor depth are given at the bottom left and right of each panel, respectively. The DART data are band-pass filtered with lower and upper cutoff periods of 4 mins and 4 hrs. The panels are sorted by increasing distance from top left to bottom right.

5.4.4.4 Additional source contributions

It is apparent from the DART signal comparison that not all of the significant peaks in the signal are predicted by either the OWC or TWC model forced by an atmospheric wave. Indeed, the pressure signal associated to the largest waves in, e.g. DART sensors 46403, 32413 32402 or 32403 are underestimated by the current modelling parameters. However, it is noted that the arrival times of these larger waves are in line with classical tsunami travel time predictions (see, e.g. Gusman & Roger Reference Gusman and Roger2022). This suggests that an initial $\mathscr {G}$-mode contribution at the volcano (which would represent, e.g. mass ejection or terrain collapse following the explosion, Matoza Reference Matoza2022) is required to achieve a quantitative prediction of the observed late-arriving signals in the DART data. Simulations carried out by injecting varying energy levels into an initial $\mathscr {G}$ mode with a Gaussian spatial support (not shown here) suggest that some of the main extrema can be recovered; this is also supported by numerical exploration by, e.g. Lynett et al. (Reference Lynett2022) although in the context of a OWC model. However, modelling the complex mechanisms governing the initial $\mathscr {G}$-mode generation exceeds the scope of this paper and will be explored in future work.

5.4.4.5 Refraction problem

As observed in the numerical simulation results, transitions from shallow to deep water (in particular at the Tonga–Kermadec trench) are zones of energy transfer from $\mathscr {A}$ modes to $\mathscr {G}$ modes. This transfer can be understood by formulating the refraction problem using the derived eigenmodes of the TWC system. First, the interaction of an $\mathscr {A}$ mode with a sharp change in bathymetry is considered. This situation arises, for example, when the sea floor rises from the deep sea to a continental shelf over a spatial scale that is small compared with the $\mathscr {A}$-mode wavelength (e.g. given a wavelength of about $10^3\ \mbox {km}$, bathymetry changes occurring over a few tens of kilometres are considered ‘sharp’). As a first approximation, the sharp depth change is considered as a discontinuity in the water depth separating two uniform, quiescent regions to the left and right of the discontinuity, respectively, with depths $H_L$ and $H_R$ (see figure 14 for scenario notations and illustration). Here $H_L>0$ and $H_R>0$ is enforced such that the step change is always submerged. To the left, a lone downstream-propagating $\mathscr {A}^+$ mode is imposed (i.e. the amplitude of the upstream $\mathscr {T}$ mode and $\mathscr {G}^+$ mode are set to zero). The interaction of the incident $\mathscr {A}$ mode with the discontinuity can generate an upstream-propagating $\mathscr {G}^-$ mode, a reflected upstream-propagating $\mathscr {A}^-$ mode, a downstream $\mathscr {T}$ mode, a downstream-propagating $\mathscr {G}^+$ mode and a transmitted downstream-propagating $\mathscr {A}^+$ mode. Thus, accounting for each of the contributions, the fluctuations to the left of the discontinuity can be expressed as

(5.6)\begin{equation} \hat{\boldsymbol{\mu}}_L = \alpha_{\mathscr{A}_L^+} \hat{\boldsymbol{\mu}}_{\mathscr{A}_L^+} + \alpha_{\mathscr{G}_L^-} \hat{\boldsymbol{\mu}}_{\mathscr{G}_L^-} + \alpha_{\mathscr{A}_L^-} \hat{\boldsymbol{\mu}}_{\mathscr{A}_L^-},\end{equation}

and the fluctuations to the right of the discontinuity can be expressed as

(5.7)\begin{equation} \hat{\boldsymbol{\mu}}_R = \alpha_{\mathscr{A}_R^+} \hat{\boldsymbol{\mu}}_{\mathscr{A}_R^+} + \alpha_{\mathscr{G}_R^+} \hat{\boldsymbol{\mu}}_{\mathscr{G}_R^+} + \alpha_{\mathscr{T}_R} \hat{\boldsymbol{\mu}}_{\mathscr{T}_R}, \end{equation}

The amplitude $\alpha _{\mathscr {A}_L^+}$ is prescribed, i.e. it is a known quantity. At the location of the discontinuity of water depth, continuity of the $\hat {\boldsymbol {\mu }}$ field is imposed, i.e. $\hat {\boldsymbol {\mu }}_L = \hat {\boldsymbol {\mu }}_R$. Thus, together, (5.6) and (5.7) form a linear system of five equations with five unknowns that are the amplitude of the reflected, transmitted and generated modes,

(5.8)\begin{equation} \begin{bmatrix} \hat{\boldsymbol{\mu}}_{\mathscr{G}_L^-}, \hat{\boldsymbol{\mu}}_{\mathscr{A}_L^-}, -\hat{\boldsymbol{\mu}}_{\mathscr{A}_R^+}, -\hat{\boldsymbol{\mu}}_{\mathscr{G}_R^+}, -\hat{\boldsymbol{\mu}}_{\mathscr{T}_R} \end{bmatrix} \begin{bmatrix} \alpha_{\mathscr{G}_L^-} \\ \alpha_{\mathscr{A}_L^-} \\ \alpha_{\mathscr{A}_R^+} \\ \alpha_{\mathscr{G}_R^+} \\ \alpha_{\mathscr{T}_R} \end{bmatrix} ={-} \alpha_{\mathscr{A}_L^+} \hat{\boldsymbol{\mu}}_{\mathscr{A}_L^+}.\end{equation}

By analogy with the derivation of (5.8), the system governing the case of a $\mathscr {G}_L^+$ mode interacting with a step results in an almost identical system where the right-hand side is replaced by $-\alpha _{\mathscr {G}_L^+} \hat {\boldsymbol {\mu }}_{\mathscr {G}_L^+}$. Once the systems are inverted (which is done here numerically with the standard atmosphere values from § 5.1 for steps ranging up to $H=12$ km in depth), two quantities of particular interest in understanding the historical maximum heights obtained by the simulations can be derived. The first is the $\mathscr {G}$-mode-associated water-height fluctuation generated as an $\mathscr {A}$ mode crosses a step, where the incident $\mathscr {A }$ is characterised by its pressure fluctuation at sea level. The second is the generated $\mathscr {G}$-mode-associated water-height fluctuation (relative to the local depth) as a result of an incident $\mathscr {G}$ mode interacting with a step, where the incident $\mathscr {G}$ is characterised by its water-height fluctuation (relative to local depth). To this end, two ratios are formed, respectively characterising each of the aforementioned scenarios, i.e.

(5.9a,b)\begin{equation} \alpha^*_1 \equiv \frac{\alpha_{\mathscr{G}^+_R} \hat{\eta}_{\mathscr{G}^+_R} }{ \alpha_{\mathscr{A}^+_L} \hat{p}_{\mathscr{A}^+_L}} \frac{{\rm \pi}'}{p'(\eta_1)}, \quad \alpha^*_2 \equiv \frac{\alpha_{\mathscr{G}^+_R} \hat{\eta}_{\mathscr{G}^+_R} }{ \alpha_{\mathscr{G}^+_L} \hat{\eta}_{\mathscr{G}^+_L}} \frac{H_L}{H_R} = \frac{\alpha_{\mathscr{G}^+_R} }{ \alpha_{\mathscr{G}^+_L}}. \end{equation}

The value of $\alpha ^*_1$ and $\alpha ^*_2$ are illustrated in figure 15 over a range of forward-facing and backward-facing steps (this orientation is decided by the direction of propagation of the incident wave). Qualitatively, the $\alpha ^*_1$ map highlights a non-symmetric behaviour when an $\mathscr {A}$ mode crosses a step: a backward-facing step generates a water-height trough whereas a forward-facing step generates a water-surface peak. Furthermore, unlike the equivalent problem for the OWC system (which contains a singularity when the set $\mathscr {A}^\pm _{owc}$ speed matches the water-layer gravitational speed; see, e.g. Vennell Reference Vennell2007), the crossing of the critical height in both step configurations happens continuously. For the incident $\mathscr {G}$-mode scenario, the $\alpha ^*_2$ map shows that the transmitted $\mathscr {G}$ is attenuated by a backward-facing step and amplified by forward-facing step.

Figure 14. Illustration of the refraction problem formulated at a forward-facing step in quiescent flow. In this scenario where no downstream propagating $\mathscr {G}_L^+$ and no upstream propagating $\mathscr {A}_R^-$ are imposed, an $\mathscr {A}_L^+$ mode interaction with a step change in depth is associated with two upstream-propagating modes ($\mathscr {A}_L^-$ and $\mathscr {G}_L^-$), one stationary $\mathscr {T}_R$ mode and two downstream-propagating modes ($\mathscr {A}_R^+$ and $\mathscr {G}_R^+$). The black arrows indicate the direction of propagation. Subscripts $L$ and $R$ correspond to left and right conditions, respectively.

Figure 15. (a) Qualitative world bathymetry map (cream to dark blue for shallow to deep water) with numbered location of sharp changes in sea floor depth. (b,c) Refraction maps for coefficients $\alpha ^*_1$ (b) and $\alpha ^*_2$ (c) over a range of depth ranges. The numbered values are the features identified in the top panel with the step orientation (i.e. forward- or backward-facing step) as seen when travelling along an isoazimuth line departing from the HTHH volcano to its antipode. The critical height $H_c$ is indicated for both axes by the dashed grey lines. For reference, the $5$ km water-depth line is indicated by dashed white lines.

To illustrate how these maps can be used to understand the aftermath of the Tonga event, bathymetry features that fit the approximations of this one-dimensional approach are identified. To do so, sharp changes along great circles (from the volcano to its antipode) that correspond to ridges that are locally approximately normal to the wave front along its direction of propagation are manually selected. To assign a ‘left’ and ‘right’ depth, the water depth either side of the manually selected features is averaged along the great circle over a distance of up to $5\times 10^2\ \mbox {km}$ (less if land is reached first). The retained features are numbered on the map in figure 15. A prominent feature of the maximum historical water-height map in figure 9 is the significant $\mathscr {G}$-mode generation as the $\mathscr {A}$ crosses the Tonga and Kermadec trenches to the southeast of the volcano. As predicted by the map in figure 15, a large amplitude depression of the sea surface is generated (due to the atmospheric fluctuations being the largest in magnitude close to the source) and propagates away to the south and east, as shown in figure 11 (row 1 and row 3). When the wave reaches locations such as the southern tip of South America and the step down off the coast of Argentina (row 4), the ground pressure fluctuation $p'(\eta _1) \sim 1.5\ \mbox {hPa}$ will, according to the refraction map, generate $\mathscr {G}$-mode associated $\eta _1' \sim -1$ cm. The generated $\mathscr {G}$ modes interact with each other, as seen in figure 11 (row 4) and in supplementary movie 4, to generate a local streak of maximum water-height disturbance.

5.4.4.6 Energy injection considerations

Estimates of the amount of either the total energy released or the energy released by the volcano into a subset of layers from empirical correlations and extrapolation from observations have yielded a wide range of values spanning multiple orders of magnitude from $10^{16}\ \mbox {J}$ to $10^{19}\ \mbox {J}$ (see, e.g. Astafyeva et al. Reference Astafyeva, Maletckii, Mikesell, Munaibari, Ravanelli, Coisson, Manta and Rolland2022; Díaz & Rigby Reference Díaz and Rigby2022; NASA 2022; Wright et al. Reference Wright2022; Yuen et al. Reference Yuen2022). The complexity of the volcanic explosion process (which included multiple smaller explosions before and after the main one, which is the focus of this study) prevents a detailed account of energy injection per layer; however, numerical studies allow for testing of a range of energy injection per layer and per mode. Unlike in the OWC case where the forcing is applied at the surface (i.e. without taking into account the thickness of the perturbed atmospheric layer), the TWC model allows for an energy balance to be performed over the volume of each layer to determine the amount of energy injected by the volcano source condition into the ocean and the atmosphere, respectively. As the source imposes a superposition of equal amplitude $\mathscr {A}^+$ and $\mathscr {A}^-$ modes, neither layer is forced with a velocity fluctuation. An accounting of the energy injection per layer is proposed as follows.

1. (a) Atmospheric-layer energy injection: the energy injected into the atmospheric layer, $E_{atm}$, can be computed by integrating the absolute mass-specific Favre-averaged internal-energy fluctuation over its spatial and temporal support. The EoS of the atmospheric layer, ${{\rm \pi} } = (\gamma -1 ) {\varrho } \tilde {e}$, and the relation between $\mathscr {A}$-mode-induced average-pressure and average-density fluctuations can be used to express the injected energy as

   (5.10)\begin{align} E_{atm} \approx \dfrac{1}{\tau} \int_0^\tau \int_S h_0 \varrho_0 | \tilde{e}' | \,\mathrm{d}S \,\mathrm{d}t \approx \dfrac{1}{\tau} \int_0^\tau \int_S \dfrac{h_0}{(\gamma -1) \varrho_0} \left( 1 - \dfrac{ {{\rm \pi}_0} }{\varrho_0 \mathscr{C}_0 ^2} \right) | {{\rm \pi}'} | \, \mathrm{d}S \, \mathrm{d}t. \end{align}
   Evaluating this integral with dimensional quantities for the conditions of the TWC simulation, detailed in § 5.4.2, yields
   (5.11)\begin{equation} E_{atm}^* = 1.7 \times 10^{16}\ \mbox{J} . \end{equation}

2. (b) Water-layer energy injection: the energy injected into the water layer can be computed by considering the potential energy contribution of the $\mathscr {A}$-mode-induced $\eta _1$ fluctuations. The integral over the source's spatial and temporal support can be expressed as

   (5.12)\begin{equation} E_{water} = \dfrac{1}{\tau} \int_0^\tau \int_S \rho_w g H | \eta_1' | \,\mathrm{d}S \,\mathrm{d}t . \end{equation}
   Evaluating this integral with dimensional quantities for the conditions of the TWC simulation, detailed in § 5.4.2, yields
   (5.13)\begin{equation} E_{water}^* = 4.0 \times 10^{14}\ \mbox{J} . \end{equation}

3. (c) Total energy injection: the total energy injected into the coupled system, $E_{tot}^* \equiv E_{atm}^* + E_{water}^*$ via $\mathscr {A}$-mode contributions is in line with the lower end of the volcanic energy release estimates. Note that this does not include any $\mathscr {G}$-mode contributions discussed above nor any seismic energy. Future work will include a study of the required $\mathscr {G}$-mode energy contribution to recover the missing peaks observed in the DART signal data.

5.5. A posteriori assumption verification

In this section a posteriori checks of the validity of the assumption made while deriving the model and linear theory are proposed.