Detailed examples
HungaTonga Hunga Ha’apai volcano explosion
On January 15th 2022, the explosion of the HungaTonga Hunga Ha’apai caused unexpectedly large wave amplitudes around the world. Following the confusion surrounding warnings issued by the Japanese Meteorological Agency (JMA) witnessed by some of the authors on Okinawa, Japan, we decided to investigate the physics involved in the tsunami generation mechanisms. Rapid first attempts were made using a oneway coupled system (Euler equations for the atmosphere, shallow water equations for the water layer; the results were used to generate the animation on the index page of this website) which led us to a more thorough analysis of the theory and the simulations of the newly developed twoway coupled model illustrated hereafter. This page details the elements that go into simulating the January 15 \(^{th}\) explosion using dNami, from setting up the symbolic equations to postprocessing the output. A grid convergence study and a comparison to available measurements are also provided. All the scripts necessary to run the case are located in the exm/2d_tonga_TWC
directory on the tonga
branch (please checkout this branch to obtain the source files). This case was developed for and included as part of a published paper titled ‘Twoway coupled longwave isentropic oceanatmosphere dynamics’ (doi:10.1017/jfm.2023.131).
Equations
Resulting from the derivation in JFM paper, the governing equations to be integrated using dNami are :
where the \(\nabla^R\) operator is expressed on a spherical shell at a constant radius \(R\), \(\eta_1\) is the water height, \(h_w\) is the water depth, \(\mathbf{V}\) is the water velocity field, \(\rho_w\) is the water density, \(g\) is acceleration due to gravity, \(\rho_a\) is the average atmospheric density, \(\Psi \equiv h_a^{1} \nabla^R ( h_w \mathbf{V} + h_a \mathbf{U})\), \(h_a\) is the atmospheric thickness, \(p_a\) is the average atmospheric pressure and \(c_a^2 \equiv \gamma p_a/\rho_a\). For more details on the derivation of the governing equations, please see [JFM paper]. As in the paper, \(\theta\) and \(\phi\) are used to denote the colatitude and longitude, respectively. In practice, to avoid spurious noise generation, the governing equations are expanded such that no numerical derivative of a trigonometric function appears (all derivatives are computed analytically beforehand). These are the form of the equations that appear in the equations.py
file. Furthermore, to lighten the syntax, a number of intermediate functions are included in the functions.py
file to generate the pseudocode for the spherical gradient and divergence terms that appear. The specification of the governing equations is then:
#  NB: iAtmos switches atmospheric terms on/off
#  Filling advectiontype terms
#  water
Adv['eta1'] = 'vt *' + grads('hw',0) + ' + vp *' + grads('hw',1)
Adv['vt'] = 'vt *' + grads('vt',0) + ' + vp *' + grads('vt',1)
Adv['vp'] = 'vt *' + grads('vp',0) + ' + vp *' + grads('vp',1)
if iAtmos:
#  air
Adv['rhoa'] = 'ut *' + grads('rhoa',0) + ' + up *' + grads('rhoa',1)
Adv['ut'] = 'ut *' + grads('ut' ,0) + ' + up *' + grads('ut' ,1)
Adv['up'] = 'ut *' + grads('up' ,0) + ' + up *' + grads('up' ,1)
Adv['pa'] = 'ut *' + grads('pa' ,0) + ' + up *' + grads('pa' ,1)
#  Filling sourcetype terms (WARNING: written on LHS!)
#  water
S['eta1'] = 'hw * ( ' + divs('vt', 'vp') + ')'
S['vt'] = grads('eta1', 0)
S['vp'] = grads('eta1', 1)
if iAtmos:
#  Psi
psi = '1.0_wp/ha * (' + divs('(ha*ut + hw*vt)', '(ha*up + hw*vp)') + ')'
##  Coupling terms
S['vt'] += '+ 1.0_wp/rhow *' + grads('rhoa*ha',0)
S['vp'] += '+ 1.0_wp/rhow *' + grads('rhoa*ha',1)
##  air
S['rhoa'] = 'rhoa *' + psi
S['ut'] = grads('eta1a', 0) + ' + 1.0_wp/(rhoa*ha) *' + grads('pa*ha',0)
S['up'] = grads('eta1a', 1) + ' + 1.0_wp/(rhoa*ha) *' + grads('pa*ha',1)
S['pa'] = 'rhoa * ca2 * ' + psi
In addition to the primitive variable arrays, a number of stored quantities are used in this case:
The various trigonometric quantities that appear in the governing equations: \(\small{sin(\theta), cos(\theta), 1/sin(\theta)}\). These values are filled once in the
compute.py
at the start of the computation.The bathymetry/topography: \(\small{\eta_0}\). This is filled by loading values at the start of the computation.
The historical maximum of the water fluctuations, \(\small{hmax}\), and the time at which it occurs, \(\small{tmax}\). These are updated at a set interval during the computation.
The temporal forcing of the RHS related to the initial/source conditions: \(\small{Svar}\) type terms. See sections below for details.
The temporal forcing of the RHS related to the nonstationary nature of the baseflow: \(\small{rhs\_var}\) type terms. They simply contain the negative of the sum of the other quantities making up the righthand side. These are computed at the very start of the computation.
Numerics
To compute the spatial derivative, a standard 5 point, 4^{th} order centered finite difference scheme is used. For numerical stability, an optimised 13 point, 8^{th} order filter is used. The code block specifying these choices can be found below and in the genRhs.py
file.
#GenerateRHS:
append_Rhs(Adv, 5, 4, rhsname, vnamesrc_Adv, update=False, rhs=rhs, stored=True)
append_Rhs(S, 5, 4, rhsname, vnamesrc_S, update=True, rhs=rhs, stored=False)
append_Rhs(C, 5, 4, rhsname, vnamesrc_C, update=True, rhs=rhs, stored=False)
append_Rhs(Fss, 5, 4, rhsname, vnamesrc_Fss, update=True, rhs=rhs, stored=False)
#GenerateFilters:
genFilter(13, 8, len(varsolved), rhs=rhs)
Note the addition of the C
and Fss
quantities to the RHS. These are the volcano source condition and the forcing terms related to the nonsteady nature of the baseflow (see following section for details).
Preprocessing and Computation
Bathymetry
Before starting the computation, global bathymetry data (obtained from GEBCO) combined with polar ice coverage information for January 2022 (data obtained from NSIDC in polar coordinates is projected to the GEBCO data grid) is interpolated to the computational grid and written to a binary file in the dNami data format. This can then be read in parallel at runtime using the read_data()
function. The misc.py
file contains a utility for the steps involved in reading the bathymetry data and interpolating it. However, the thermodynamic fields constructed this way no longer respects the spherical boundary conditions and imposing them can lead to numerical instability. To remedy this, the fields are smoothed normal to each boundary (using a half sinwave) to the average value either side of the boundary over a thin layer near the boundary (chosen to be 10 grid points).
Verticallyaveraged atmospheric conditions
Similarly to the bathymetry data, prior to starting the computation, the underlying fields of velocity and thermodynamic quantities are interpolated to the computational grid and written to dNami binary data format. They are derived from ERA5 data for January 15th 2022 with can be found on the Copernicus Climate Data Store plateform . However, just like the thermodynamic fields, the velocity field constructed this way no longer respects the spherical boundary conditions and imposing them can lead to numerical instability. The same smoothing strategy as above is employed to remedy this.
Nondimensionalisation values
A Python class containing the values used to nondimensionalise the input quantities is provided in values.py
. Their physical significance is given in the table below.
Variable name 
Physical meaning 


Average water depth on Earth 

Characteristic gravitational wave speed (in water) 

Characteristic time scale based on 

Reference water density 

Reference ‘dynamic’ pressure based on 
Imposing spherical boundary conditions
By default, a double periodic domain in dNami imposes Cartesian geometry. Fig. 11 illustrates a configuration where the domain is shared between 8 processors (4 processors in the horizontal direction and 2 in the vertical direction). To impose spherical boundary conditions, the vertical neighbours for processors along the poles must be redefined. In addition, the horizontal swap operations are maintained whereas the vertical operations involve matching processors along the poles of the spherical domain (i.e. North or South) and flipping (both horizontally and vertically) the contents of the swapped information. These changes are illustrated for processor 3 where the bottom vertical neighbour of processor 3 is no longer 7 but 0 and the flipping of the information in the bottom halos is illustrated by the ‘L’ symbol (as seen by procs 0 and 3).
To be able to match processors onetoone along the poles, an even split in number of processors is imposed. The reorganisation of the processors (from the initial Cartesian geometry) is done in the following code block found in the compute.py
file which sets flags in the dMpi
class used in the modified dnami_mpi.py
file distributed with this case:
#  Require an even split
if with_proc[1] % 2 != 0:
if dMpi.ioproc: print('[ERROR] ydirection split must be even ...')
exit()
#  Define a new position/negative direction
pid = dMpi.procid
mod = np.mod(pid, with_proc[1])
if mod <= int(with_proc[1]/2)  1 :
#Copy
nxm = 1*dMpi.neighxm
nxp = 1*dMpi.neighxp
#Swap
dMpi.neighxm = nxp
dMpi.neighxp = nxm
dMpi.flipped = True
#Flip ud and lr
if dMpi.ibeg1 == 0 or dMpi.iend == grid['nxgb']:
dMpi.reverse = True
#  South pole procs
if dMpi.ibeg1 == 0:
mid = int(with_proc[1]/2)  1
if pid <= mid:
dp = abs(midpid)
dMpi.neighxp = mid + dp + 1
dMpi.edge = True
else:
dp = abs(mid+1pid)
dMpi.neighxm = mid  dp
#Flip ud and lr
dMpi.reverse = True
#  North pole procs
if dMpi.iend == grid['nxgb']:
mid = with_proc[0]*with_proc[1] 1  int(with_proc[1]/2)
if pid <= mid:
dp = abs(midpid)
dMpi.neighxm = mid + dp + 1
else:
dp = abs(mid+1pid)
dMpi.neighxp = mid  dp
#Flip ud and lr
dMpi.reverse = True
dMpi.edge = True
##  Velocity flipping  lat component must be opposite (in this case x component)
dMpi.sym = dMpi.reverse and dMpi.edge
if iAtmos:
dMpi.velidx = [1,4] #vt, ut
else:
dMpi.velidx = [1] #vt
Volcano source conditions
Three elements are important in imposing the volcano source:
the relative contribution of the source to each primitive variables
the spatial support
the temporal support
For the first point, the amplitude of the source for each of the primitive variables is given by the components of the eigenvectors (see JFM for details). Given the nature of the atmospheric wave, a superposition of an \(A^+\) and an \(A^\) mode is prescribed at the origin such that only the atmospheric pressure and density and waterlevel fluctuations are disturbed. The amplitudes of each of these contributions are computed in the following code block:
# Compute the relevant eigenmode amplitudes
etahat = phi*H0
rhohat = rhoa[dom] *(lam2/C02 + Cw2/C02*(beta1.))
phat = rhoa[dom]*C02*(lam2/C02 + Cw2/C02*(beta1.))
For the second point, to avoid numerical (and physical) issues with a single point source (e.g. the addition of nonphysical frequency due to the discretisation operation), the source is given a spatial support in the form of a 2D Gaussian. To ensure the Gaussian is circular (in easting/northing coordinates and not in latitude/longitude coordinates), the following code block leveraging the haversine python package is used to generate a Gaussian support based on the distance of each point to the volcano (assuming a spherical Earth).
def gauss(x0,y0,amp,rad0):
from haversine import haversine_vector, Unit
# Generate a Gaussian support based on the distance from each point to the volcano
rad0m = rad0*Re
x0d = x0*180./np.pi  90.
y0d = y0*180./np.pi  180.
xld = xloc[:,np.newaxis]*180./np.pi  90.
yld = yloc[np.newaxis,:]*180./np.pi  180.
XX,YY = np.meshgrid(xld, yld, indexing='ij')
darray = haversine_vector(list(zip(x0d*np.ones_like(XX).ravel(),y0d*np.ones_like(YY).ravel())),list(zip(XX.ravel(), YY.ravel())), Unit.METERS).reshape((nx,ny))
return amp*np.exp( (darray**2)/(2.*rad0m**2))
For the last point, the explosion is modelled as a local source with a given temporal shape (amplitude and period). The \(Svar\) terms are updated in time and used to force the RHS to give the source its temporal support. Given the shape of the observed wave and to ensure a continuous forcing (at least \(\mathcal{C}^1\)), a 5 \(^{th}\) order polynomial is constructed with a temporal support of duration \(\tau\) satisfying the following constraints: the slope and value of the function at \(t=0\) and \(t=\tau\) must be zero and the function must reach \(p^+\) at \(t=\tau/4\) and \(p^\) at \(t=3\tau/4\). Note that this term is applied as a forcing to the RHS by premultiplying the amplitude of the spatial Gaussian support (see following animation for illustration).
def update_ft(t):
# Default value
ft = 0.
a = (512.*pminus  512.*pplus)/(9.*tau**5)
b = (1152.*pminus + 1408.*pplus)/(9.*tau**4)
c = (768.*pminus  1280.*pplus)/(9.*tau**3)
d = (128.*pminus + 384.*pplus)/(9.*tau**2)
#  Force if time less than tau
if t <= tau:
ft= a*t**5 + b*t**4 + c*t**3 + d*t**2
#  Update coefficient used to force RHS
fltparam[5] = ft
dtree['eqns']['coeff'][0][1] = ft
return
If this spatiotemporal forcing is integrated in time in 1D Cartesian coordinates with a spatial scale \(\sigma_x = 50\) km and temporal scale \(\tau = 45\) min, the field of density, waterlevel and pressure fluctuations are obtained as shown in the animation below. The parameters governing the forcing allow for the wavelength and amplitude of the crest and trough of the resulting atmospheric wave to be adjusted in accordance with the observations.
Extracting station information
To compare the pressure or waterwave height at various measurement stations around the world to the result of the computation, two possibilities are available to the user: either extract the desired value at the station coordinates during the computation (coprocessing) or output restart files and extract the value from the file after the computation (postprocessing). The first approach can slow down the computation as it will unbalance the load between processors whereas the second can result in large storage requirements (i.e. restarts need to be stored at the correct frequency to reconstruct the desired temporal signal).
To implement the first option, prior to starting the time advancement, the list of stations around the world is read and each processor retains the coordinates in its subdomain as shown in the following code block.
#  Load station coordinates (each procs takes its own coordinates)
coords_atm = {}
#  Ground pressure stations
for station in glob.glob('stations/**/*.bin', recursive=True):
key = station.split('/')[1]
with open(station, 'rb') as fh:
lat = np.asarray(np.load(fh, allow_pickle=True)).item()
lon = np.asarray(np.load(fh, allow_pickle=True)).item()
#  Add coord to dic if point in proc
xp = (lat+90)/180*np.pi
yp = (lon+180)/180*np.pi
if (xp > xloc[0]0.5*dx and xp<xloc[1]+0.5*dx ) and (yp >yloc[0]0.5*dy and yp < yloc[1]+0.5*dy):
coords_atm[key] = [xp, yp]
During the time loop, the value of the relevant field is interpolated by the corresponding processor using scipy’s map_coordinates()
function (here with alias MC
) and append to the relevant file to update the time history:
for key in coords_atm.keys():
# Get coordinates from dic
xp, yp = coords_atm[key]
# Global coordinate
nxp = xp*grid['nxgb']/Lx
nyp = yp*grid['nygb']/Ly
# Proc. local coordinate
nxpl = nxp  (dMpi.ibeg 1)
nypl = nyp  (dMpi.jbeg 1)
# Interpolation
out = MC(pa[dom], [[nxpl], [nypl]], order=1)
# Write out
with open(outf+'/stations/'+key[:4]+ '_' + str(n).zfill(6) + '.bin', 'wb') as fh:
np.save(fh,np.array([ti, out.item()]))
fh.close()
RHS forcing
Due to the nonstationnary nature of the initial conditions (the underlying fields of pressure, density and velocity are not a steadystate solution to the governing equations), the RHS (without yet having introduced any perturbation related to the explosion) is computed and saved in rhs_var
variables (where var
is replaced by the corresponding component of the primitive variable array). As per the equations.py
file, this value is subtracted from the rest of the RHS at each timestep. The code block below gives these steps and also shows that these forcing terms are written to file (in the event of a restart of the simulation).
#  Get RHS forcing terms
if iGetRhs and not iRestart:
#  Swap q fields
dMpi.swap(q,hlo,dtree)
#  Compute forcing terms  ONLY DO THIS ONCE
if 'qstored' in dtree['eqns']['qvec']['views'].keys():
dn.dnamiF.stored(intparam, fltparam, data)
dMpi.swap(qstored,hlo,dtree,spherical=False)
#  Write them out if used for restarting
rhsnames = []
for key in fsskeys:
rhsnames.append( 'rhs_' + key )
if dMpi.ioproc:
print(' > Saving ...', rhsnames)
dn.dnami_io.write_data(rhsnames,0,0,dtree,fpath=outf+'force/',fname='rhs')
# ...
I/O and timeloop updates
During the time loop, at intervals governed by the value of tdup
, after the solution has been updated by the RK steps, the historical waterlevel fluctuation and the occurance time are updated by:
#  Update historical hmax and tmax
if tdmax > tdup :
idx = np.nonzero(np.abs(eta1[dom]) > (hmax[dom]+1e8) ) # add epsilon to remove noise (e.g. filtering) from time history
hmax[dom][idx] = np.abs(eta1[dom][idx])
tmax[dom][idx] = ti
tdmax = 0.
Following this, at intervals governed by the value of tdcstep
, a restart file is written and the updated values of hmax
and tmax
are output.
if tdc > tdcstep :
#  Write restart
dn.dnami_io.write_restart(n,ti,0,dtree,fpath=resf)
#  Write hmax/tmax
dn.dnami_io.write_data(['hmax', 'tmax'],n,ti,dtree,fpath=outf+'historical/',fname='hmax_tmax')
#  Reset time
tdc = 0.
Output and postprocessing
The plot_water_air.py
script provides a basic utility to plot the simultaneous water and ground pressure level fluctuations. The script reads in the list of restart files and it will output images to generate animations such as the one below.
Convergence and comparison to measurements
To demonstrate the convergence of the results with grid refinement and to illustrate how the model compares to measurements on the day, a few comparison are provided below. On the issue of grid convergence, the spatial scale of the volcano source is scale to ensure it is always ‘sufficiently resolved’ at lower grid sizes to avoid the discretisation of a pointsource or underresolved perturbation. Here, ‘sufficiently resolved’ is defined as \(\sigma_x/\Delta \geq 5\) where \(\Delta\) is the spatial resolution at the source. Therefore, the spatial scale is to be \(\sigma_x = max(\sigma_x^t, 5\Delta)\) where \(\sigma_x^t = 50\) km is the target spatial scale. This corresponds to onescale separation away from the shock propagation scale (which is also in agreement with the vertical span of the perturbation). As such, the final values of waterlevel and atmospheric pressure are only expected to converge after the target spatial scale is ‘sufficiently resolved’. The crossing of this threshold as a function of grid size is illustrated in Fig. 14 as well as the grid sizes presented below.
Historical maximum
The historical maximum of the absolute sea surface fluctuations are given at various grid sizes as an illustration of convergence of the results with mesh refinement. Note however that the highest resolution available for the bathymetry grid is 86400x43200 thefore smaller scale features will continue to appear in the final solution until that grid resolution is reached.
Pressure measurements
Four sources are used to gather pressure measurements from around the world: Weather News Inc. from Japan (kindly provided a dataset upon request, see here), Iowa State University hosts 1 minute data from ASOS sites in the USA (IASTATE), Sensor Community archive data (SensCom) and a signal from the University of Reading (UoR). The location of each of these sensors is given in Fig. 16
Values of ground level atmospheric pressure fluctuation at the corresponding locations are obtained from dNami by bilinear interpolation of the pressure on the computational grid (either at runtime or afterthefact). The result is shown in (where a tag of the same colour as in Fig. 16 has been added to identify to which group the measured signal belongs).
To demonstrate convergence of the pressure signal with grid refinement (and corresponding convergence of the signal’s spatial support), results from computations on increasingly refined grids are presented in Fig. 18.
DART tide data
Data from Deepocean Assessment and Report of Tsunamis (DART) probes was gathered from the NOAA. The location of the probes is shown in Fig. 19. In deep water, where the hypotheses involved in deriving the shallow water equations are satisfied, a quantitative comparison can be made between measurements and results from the theory. It should be noted however that the water height fluctuation given by the DART probes does not take into account atmospheric pressure fluctuations.
Similarly to what is done for the pressure sensors, signals from the computation are extracted at each DART probe location and compared to the realworld measurements in Fig. 20. For comparison, the DART data is filtered with a bandpass filter with cutoff periods of 4 minutes and 4 hours (to remove wave noise and tide related content). In addition, for a compatible comparison with what is provided by DART, the ocean bottom pressure fluctuations are compared (rather than the waterheight fluctuations). NOAAsourced DART data is converted back to bottom pressure fluctuations using the constant 670mm/psia given by DART documentation.
To demonstrate convergence of the waterlevel fluctuation with grid refinement, results from computations on increasingly refined grids are presented in Fig. 21. Note that smaller wavelength contributions will continue to appear with grid refinement until the maximum resolution of the bathymetry is reached (roughly 400m/direction).
Waterlayer only
To illustrate the impact of the coupling and contribution of the atmospheric wave, a case without the atmospheric layer is run. Only the water layer is advanced in time with a source condition which uses only the \(\eta_1\) contribution detailed above. The resulting historical maximum waterheight fluctuation is given in Fig. 22 over the same 18 hour period. Note how, given the slower speed of gravity wave propagation in the water and the necessity to circumvent land, the perturbation does not travel as far, barely making it past South America and does not reach Africa. Furthermore, the energy injected into the water layer, without the atmospheric coupling, is not enough to explain the observed global tsunamis in the wake of the volcano eruption.