Adjoint Synthesis for Trans‐Oceanic Tsunami Waveforms and Simultaneous Inversion of Fault Geometry and Slip Distribution

Tsunamis propagate over long distances and can cause widespread damage even after crossing ocean basins. Prediction of tsunamis in distant areas based on observations near their sources is critical to mitigating damage. In recent years, the accuracy of numerical models of trans‐oceanic tsunami propagation has improved significantly due to the incorporation of effects such as the solid earth response to tsunami loading and wave dispersion. However, these models are computationally expensive and have not been fully utilized for real‐time prediction. Here, we derive the adjoint operator for the linear set of equations describing deep‐ocean tsunami propagation and show how a pre‐computed database of adjoint states can achieve rapid synthesis of tsunami waveforms at target sites from nonpoint arbitrary tsunami sources. The adjoint synthesis method allows for an exhaustive parameter search for tsunami source estimation. A method for simultaneous inversion of fault geometry and slip distribution using adjoint synthesis with Sequential Monte Carlo method was proposed and applied to the 2012 Haida Gwaii earthquake tsunami. The influence of model accuracy and the amount of observed data on the estimation of tsunami sources and waveforms was examined. It was found that with a highly accurate propagation model, using only a limited amount of observed data produced source and waveform estimates very similar to the final models obtained with much larger data sets. The final inferred fault model involved megathrust slip distributed between the Haida Gwaii trench and the Queen Charlotte fault. The proposed method can also quantify the uncertainty of the waveform forecasts.


Introduction
Tsunamis are generated mainly by water surface disturbances caused by large-scale seismic motions, and their source areas range from tens to thousands of kilometers.Tsunamis transmit energy to great distance and cause damage over wide areas.The 1946 Aleutian Islands earthquake tsunami struck Hilo, Hawaii, approximately 3,700 km away from the source area, about five hours after the quake, killing 173 people.This event led to the recognition of the importance of early warnings and the establishment of the Pacific Tsunami Warning Center in 1949.Subsequently, trans-oceanic tsunami disasters have continued, such as 2004 Indian Ocean, 2010 Chile and 2011 Japan tsunamis.
Tsunami waveforms have long been used to better understand how these destructive waves are generated.Inversions of coastal tide gauge data, typically using the first cycle of the direct tsunami arrival, yielded important

10.1029/2024JB028750
Key Points: • An exact reciprocity approach for linear deep ocean tsunami waveforms is established by defining a tsunami adjoint operator • A method for simultaneous estimation of fault geometry and slip distribution using a sequential Monte Carlo method is proposed • Exhaustive parameter search enabled by adjoint synthesis improves waveform prediction accuracy for trans-oceanic tsunami Supporting Information: Supporting Information may be found in the online version of this article.information on how large megathrust earthquakes of the twentieth century generated large tsunamis (e.g., Satake, 1987Satake, , 1989Satake, , 1993;;Satake & Kanamori, 1991).Precision of earthquake rupture models determined with tsunami waveform inversion improved with the development and deployment of ocean-bottom pressure gauges in the late 1990s and early 2000s (Aoi et al., 2020;Bernard & Titov, 2015), since direct arrivals of deep ocean tsunamis at regional distances appeared to be accurately modeled using linear long-wave theory (Baba et al., 2009;Fujii et al., 2011;Lay et al., 2005Lay et al., , 2014;;Satake et al., 2013;Yamazaki et al., 2018;Yue et al., 2014Yue et al., , 2015)).However, it was generally thought that far-field stations provide limited improvements to the tsunami source.
The 2010 Chile and 2011 Japan tsunamis were the first large events that were well-recorded by ocean-bottom pressure gauges in the far field, and studies comparing the observed and predicted far-field waveforms noted systematic differences (Bai et al., 2015;Fujii & Satake, 2013;Grilli et al., 2013), including arrival-time discrepancies of up to 15 min in the far field.Tsai et al. (2013) and Watada (2013) showed that these discrepancies could be explained by considering the deformation of the solid Earth and seawater compressibility and density stratification, which we describe below as secondary effects.Subsequently, Baba et al. (2017) showed that these effects, as well as Boussinesq dispersion, can all be significant in the far field.Baba et al. (2017) showed that accounting for these effects, which we describe in more detail below, can yield excellent fits to tsunami waveforms of several hour's duration observed in both near and far field.
In this paper, we explore how accurate simulation of transoceanic tsunamis can allow the information contained in long (1-6 hr) tsunami wavetrains to be exploited for more accurate estimates of tsunami sources and forecasts of observed waveforms.In particular, we consider how accurate source characterization from near-field, extended tsunami wavetrains (e.g., 1 hr duration) can be used to accurately simulate far-field tsunami waveforms, and how these improvements in accuracy of both tsunami sources and modeled far-field tsunamis might be achieved in real time, by analyzing the progressive improvement in accuracy that can be achieved at various post-event observation times.
Improvements in the accuracy of real-time estimates of tsunami sources and waveforms are potentially important for tsunami warning systems, but such systems also require that computations and analysis are robust and efficient (e.g., Gica et al., 2008;Tsushima & Ohta, 2014).While forward tsunami computation that includes the secondary effects and Boussinesq dispersion mentioned above can be very computationally demanding (Baba et al., 2017), we develop an adjoint operator for linear tsunami propagation that allows reciprocity to be exploited in tsunami source inversions, and allows the construction of an adjoint state database for tsunami waveforms that can facilitate very efficient computation of tsunami waveforms from a designated source area, like a subduction zone, to a limited number of deep-ocean, sea level observation sites.Because this adjoint state database is developed for sea level displacements in the source area, it can be used for any tsunami source, include earthquakes on faults of arbitrary geometry or submarine landslides.
Adjoint operators have been used for tsunami full waveform inversion.The initial water level of the tsunami source was estimated by defining a misfit function between the observed and estimated waveforms and minimizing it through the gradient descent (Maeda, 2023;Pires & Miranda, 2001) or conjugate gradient method (Xie et al., 2023;Zhou et al., 2019).The adjoint operator was used to efficiently compute that gradient.In these studies, adjoint equations are derived for shallow water equations in a Cartesian coordinate system that do not take into account the secondary effects and Boussinesq dispersion.
In this study, the adjoint equations are derived for the shallow water equations in spherical coordinates that take into account the secondary effects and Boussinesq dispersion, which are expected to reproduce trans-oceanic tsunami propagation with higher accuracy.Adjoint operators are used to create an adjoint state database, which is equivalent to a densely gridded Green's function database, rather than to compute the gradients of misfit functions.Since analysis of accuracy in both estimated sources and forecasted waveforms requires robust quantification of uncertainty, we have also developed a sequential Monte Carlo (SMC), or particle filter (e.g., Buchholz et al., 2021;Nguyen et al., 2016), approach to tsunami source inversion that utilizes the adjoint state database to produce an ensemble of tsunami source estimates that can be used for a probabilistic description of earthquake-generated tsunamis having spatially varying slip on faults of arbitrary orientation.We demonstrate how this new method can be used with an adjoint state database developed for the 2012 Haida Gwaii earthquake (USGS, 2012).We show how this approach estimates a source model that achieves a better fit to observed tsunami waveform data than has been obtained in previous studies, that it provides not only an optimal source model but also an assessment of uncertainty in the source and tsunami waveforms predicted using it, and that its efficient use of the adjoint database makes it sufficiently rapid for the new approach to be considered for use in tsunami warning systems.

Accurate Ocean Basin Tsunami Propagation and Adjoint State Source Inversion
In this section we describe early and recent developments in the theory of deep ocean tsunami propagation that allow transoceanic tsunami to be modeled with high accuracy, and the development of Bayesian approaches to source inversion for earthquake-generated tsunami that account for nonlinear dependence on fault parameters and robust assessment of model uncertainty.

Gravitational Loading, Seawater Compressibility, and Dispersion
Although tsunami propagation has typically been treated as a surface gravity wave of an incompressible fluid over a rigid bottom, highly accurate observational data provided an opportunity to examine secondary effects, such as fluid compressibility, bottom deformation, and gravity potential change.In a pioneering theoretical study, Matsuzawa (1950) formulated surface gravity waves of a compressible fluid layer on a homogeneous, elastic halfspace and obtained a dispersion relation equation.Nakamura (1961) calculated dispersion relations using typical crustal parameters for a water depth of 4,500 m and showed that the compressibility of seawater and elastic deformation of the substratum were taken into account to reduce the maximum phase and group velocities of a tsunami by 1%.Tsai et al. (2013) showed that the difference between the observed tsunami arrival times and the predictions of conventional tsunami models can be quantitatively explained by the effects of the density stratification of seawater and elastic deformation of the solid earth based on a realistic multi-layered earth model.Yue et al. (2014) proposed a method to correct the phase delay of the tsunami waveforms at observation points by considering these two effects and non-shallow water dispersion (Mei, 1989).The amount of phase correction is a function of water depth.To obtain the overall phase correction amount on an actual bathymetry, they approximated the tsunami propagation path with a great circle, divided the path into many segments, and summed the phase deviations estimated from the depths of each segment.Watada et al. (2014) showed that, by additional consideration of the effect of geopotential variations associated with tsunami motion, the propagation characteristics of observed tsunamis at different frequencies can be better reproduced.They propose another phase correction method.The method assumes that the horizontal scale of the ocean depth variation is much larger than the tsunami wavelength and calculates phase corrections for a reference ocean depth and the raypath length from the source to an observation point.The raypath length is approximated by the great-circle length in Watada et al. (2014).Ho et al. (2017) proposed an efficient way to determine the raypath and its length on the actual bathymetry and used them for the phase correction.The raypath is determined as the path that minimizes the sum of the tsunami propagation time from a source and the time from an observation point.These phase corrections are applied to the waveforms simulated by a conventional long wave model.It improves the accuracy of waveform reproduction without increasing the computational cost.Both of these two correction methods are based on the raypath and cannot take into account the secondary effects, such as seawater density stratification, elastic earth deformation, geopotential variation and non-shallow water dispersion, on two-dimensional tsunami propagation including ray convergence and divergence, multiple reflections and diffraction over non-uniform bathymetry.Allgeyer and Cummins (2014) developed a numerical simulation model to reproduce the effects of elastic loading and seawater density stratification in actual bathymetry.This model allows for the consideration of these secondary effects on two-dimensional tsunami propagation over non-uniform bathymetry.Baba et al. (2017) improved this model by further incorporation of the effects of Boussinesq dispersion and gravitational potential variation, and the waveforms of distant tsunamis were reproduced very accurately using fine resolution modeling.
The correction methods of phase delay over the raypath have been applied to tsunami source inversions and waveform hindcasts (e.g., Fujii et al., 2021;Gusman et al., 2015;Ho et al., 2017;Lay et al., 2014;Yue et al., 2014).Yue et al. (2014) applied the correction method to the 2010 Maule, Chile earthquake and tsunami and revealed localized fault slip near the trench.They concluded that such correction is essential for source inversion using far-field tsunami waveforms.However, the numerical model proposed by Baba et al. (2017) that can consider those three effects and Boussinesq dispersion without path approximation has as yet not been applied to tsunami waveform inversion, even though Boussinesq dispersion has been shown to be important for precise inversion of near-field tsunamis (e.g., Gusman, Mulia, et al., 2016;Hossen, Cummins, Dettmer, & Baba, 2015;Li et al., 2016;Saito et al., 2010), probably because the computational cost is too high.

Tsunami Adjoint State Database and Finite Source Inversion
Tsunami source inversion is generally a nonlinear optimization problem that requires many repeated simulations to evaluate how the waveforms depend on source parameters like fault orientation and size.This is a major challenge when computational cost is high, such as when secondary effects and/or Boussinesq dispersion is considered, so various approaches are often used to linearize the problem.These often use a database of precomputed, unit-source tsunamis, or "Green's functions," that facilitate efficient modeling of tsunami waveforms.
One of the most popular linear parameterizations of the tsunami source is a finite fault approximation, where slip is distributed over a rectangular, planar surface (e.g., Fujii et al., 2011;Lay et al., 2005Lay et al., , 2014;;Satake, 1987Satake, , 1989Satake, , 1993;;Satake & Kanamori, 1991;Satake et al., 2013;Yamazaki et al., 2018;Yue et al., 2014Yue et al., , 2015)).The parameters such as position, strike, dip, and spatial extent, which are nonlinearly related to tsunami waveforms, are given in advance, usually based on seismic inversion analysis with empirical scaling laws and information from geological surveys and earthquake catalogs.The amount of slip, which has a linear relationship with tsunami waveforms, is used as a source model parameter to be estimated.The fault is divided into small subfaults, and the tsunami waveforms at the observation points are evaluated by tsunami simulation using unit slip in the strike and dip directions in each subfault.These tsunami waveforms are regarded as impulse responses, or Green's functions.Observed waveforms can be modeled by a linear combination of these Green's functions, with linear coefficients that correspond to the amount of slip in each direction on each segment.This linearization is so effective in reducing computational cost that it has been implemented in a part of NOAA's Short-term Inundation Forecast for Tsunamis (SIFT) (Gica et al., 2008), which uses a database of such Green's functions calculated for circum-Pacific subduction zones.However, there are some limitations due to the use of a priori information that is not varied during the inversion, such as subfault discretization.Nakata et al. (2019) showed that waveform predictions can be sensitive to horizontal position which has a nonlinear relationship with observed waveforms, and that nonlinear optimization through a grid search can greatly improve the accuracy of tsunami waveform prediction.
A different type of linearization is the direct parameterization of sea surface displacement as the tsunami source.Unit tsunami sources with unit height and simple shape, such as a truncated pyramid (Tsushima, et al., 2009), Gaussian (Tsushima et al., 2012) or raised cosine (Dettmer et al., 2016) are equally spaced, and their linear combination is used to model the initial water surface disturbance as a tsunami source.This method has advantages and disadvantages similar to the inversion of seismic waveforms for finite fault models of earthquake rupture (e.g., Hartzell & Heaton, 1983;Ide, 2007).Because the subfault grid may contain hundreds of points, such studies often take advantage of the reciprocity theorem for the elastodynamic wave equation, which under certain conditions ensures the equivalence of simulations that exchange source and station positions (e.g., Aki & Richards, 2002).This greatly reduces the number of computations required for source inversions because it allows the number of simulations required to be made proportional to a relatively small number of stations, rather than the much larger number of source grid points (e.g., Eisner & Clayton, 2001;Furumura & Maeda, 2021;Graves & Wald, 2001;Zhao et al., 2006).
A similar approach has been taken in tsunami source studies (Hossen, Cummins, Roberts, & Allgeyer, 2015;Mulia et al., 2018) using an approximate reciprocity principle developed for the linear shallow-water equations by Korolev (2011).This approximate form of reciprocity is valid for integrals of sea surface displacement over finite areas containing source and station locations, respectively, under the assumption that the water depth inside each area is approximately constant (Korolev, 2011).In practice, this condition is not always met, because tsunami sources range over scales of 10-1,000 km, and can occur in an area with complex bathymetry, such as an ocean trench.For example, for the 2012 Haida Gwaii earthquake and tsunami shown in later chapters, most of the fault slip is thought to be concentrated in an area up to 30-40 km from the Haida Gwaii trench (e.g., Gusman, Sheehan, et al., 2016;Lay et al., 2013), where water depths vary from 1,000 to 3,000 m.
Green's function reciprocity is exact for any set of linear equations that are self-adjoint (Stakgold, 1979), such as linear elastodynamic equations.This cannot be ensured for the linear shallow-water equations since they are not self-adjoint.However, as we show below, exactly the same computational efficiency that is achieved in the seismic case using point-source reciprocity can be achieved for tsunami simulations by using the adjoint operator of the system of equations for linear deep-ocean tsunami propagation.
This study proposes an adjoint operator method that does not rely on the approximations of constant water depth in the source and observation areas and allows the evaluation of perfectly matched back-propagated waveforms from nonpoint sources in any bathymetric field.Any extended source can then be synthesized by a linear superposition over a collection of such point sources.An accurate backward waveform evaluation method is possible by deriving adjoint model equations corresponding to the usual forward model equations.Here, we derive the adjoint equation from the tsunami equations with the effects of seawater density stratification, elastic loading, geopotential variation, and Boussinesq dispersion.By performing adjoint simulations for each observation point and storing the time series of the adjoint state as a pre-computed database, it is possible to efficiently evaluate the waveform at the observation points from arbitrary nonpoint sources over a source area with arbitrary variations in bathymetry.The number of simulations required to create the database can be reduced from number of computational grid cells in a potential tsunami source area to the number of observation sites.For example, Tsushima et al. (2009) created a Green's function database with 540 arc s grid spacing and performed tsunami source inversion off the Tohoku coast using observed waveform data from 18 sites.The database was constructed by running 1,698 forward simulations with a maximum computational grid size of 20 arc s.Adjoint simulations of 18 runs can generate a Green's function databases with 20 arc s grid spacing in this case.If the same densely grided Green's function database were created by forward simulations, about 1.2 million (∼1,698 × (540/20) 2 ) simulations would be required.Thus, adjoint simulation is advantageous for creating densely gridded Green's function databases.Tsunami waveform synthesis using precomputed databases of Green's functions of unit source sea surface displacements has been used in previous studies (Dettmer et al., 2016;Hossen, Cummins, Dettmer, & Baba, 2015;Tsushima et al., 2009).These are similar to the SIFT database described above (Gica et al., 2008), except that the SIFT database is based on predefined fault segments, with geometries conforming to subduction zone megathrusts.The method using adjoint state, point-source sea surface displacement Green's functions can synthesize accurate, deep ocean tsunami waveforms corresponding to any tsunami source represented at the resolution of the computational grid used in the database calculation.
In this paper, we propose a method to simultaneously estimate the geometry and slip distribution of a tsunamiproducing fault using the Adjoint Synthesis and SMC.Simultaneous estimation is expected to avoid the accuracy loss caused by predefinition of fault geometry.The proposed method is evaluated using the 2012 Haida Gwaii earthquake tsunami as a case study with respect to its performance in predicting the waveforms of distant tsunamis.

Adjoint Waveform Synthetics for Distant Tsunami
We begin with the tsunami model proposed by Baba et al. (2017).The governing equations of continuity and motion for a non-linear dispersive wave model including the effects of seawater density stratification, elastic loading, and gravitational potential change, are expressed as follows: (2) Journal of Geophysical Research: Solid Earth 10.1029/2024JB028750 TAKAGAWA ET AL.
where η = η 0 + ζ, η 0 is the difference in sea level from its value at rest, ζ is the displacement of seafloor due to tsunami loading relative to the geoid, M and N are the depth-integrated fluxes along the longitude φ and colatitude θ directions, t is time, R is the Earth's radius, h is the depth of the ocean at rest, ρ h and ρ avg are the sea water density at the seafloor and the average along the vertical profile, respectively, g is the gravitational acceleration, n is Manning's roughness coefficient, and f is the Coriolis parameter that is expressed as follows.
where Ω is the rotation rate of the earth.The second and third terms on the left-hand sides of Equations 2 and 3 are advection terms.The third terms on the right-hand sides represent the bottom friction.It is well known that these nonlinear terms have a significant impact on propagation in shallow waters, but their impact is small in deep oceans.Governing equations that take into account higher-order dispersion terms have also been proposed (e.g., Lynett, 2006).If the tsunami wavelength is much longer than the horizontal scale of the variations in bathymetry, the spatial derivative of the sea depth appearing in the high-order dispersion terms can be neglected and Equations 2 and 3 are obtained (cf., Saito, 2019).Baba et al. (2017) simulated the 2011 East Japan earthquake tsunami and showed the effect of each of the nonlinear terms on the waveforms recorded at DART buoys in deep oceans.
Comparing the results of linear dispersive wave equations and nonlinear dispersive wave equations, there is no apparent difference in the waveforms in the first hour from the tsunami arrival time, although waveform differences after this time can be seen due to the effect of later-arriving tsunami energy that has passed through shallow waters.This difference will be discussed in the case study section.Neglecting these nonlinear terms and assuming that the water level fluctuation is sufficiently small relative to the water depth, the following linear dispersion equations (e.g., Saito, 2019) are obtained.
where Φ is a variable introduced by Shigihara and Fujima (2014) to efficiently calculate the dispersion terms, and is defined as follows: In the evaluation of ρ h , ρ avg , and ζ, we follow Allgeyer and Cummins (2014) and Baba et al. (2017).
Equations 5-7 can be simplified using a vector x = (η, M, N) ⊤ and a linear operator L.
Here, we define an inner product as follows: (φ, θ), N src (φ, θ) at time zero, the water level at an observation point (φ obs , θ obs ) at a given time t is expressed as follows: where and δ is the Dirac delta function.The source function x src corresponds to an instantaneous displacement of sea surface and change of flow fields.The observation function x obs is an operator to retrieve the value of the water level at longitude φ obs , colatitude θ obs , and time t.The fact that the second and third rows of x obs are zero does not mean that these components are assumed to be zero, but simply means that only the water level component of Lx src is used in the calculation of Equation 11, which can be rewritten as: where The Riesz representation theorem (Stakgold, 1979) states that the adjoint of L, L † , is uniquely defined and satisfies the following relationship for any x and x′: The concrete expression of the adjoint operator is obtained by taking the partial integral, as shown in Text S1 of Supporting Information S1.The derived adjoint equations are as follows: Journal of Geophysical Research: Solid Earth 10.1029/2024JB028750 TAKAGAWA ET AL.
where the primed variables denote adjoint parameters.These equations are consistent with Equations 5-7 under the following variable transformations.
This consistency holds while neglecting the spatial derivative of the sea depth in the dispersion term.This approximation is also used as described above in the derivation of dispersive wave equations, such as Equations 2 and 3 (cf.Eq. 6.78 in Saito, 2019).Denoting this variable transformation of Equation 20 by P, Equation 11 can be calculated as follows: where LPx obs (0)dt ( 22) And P 1 indicates the inverse variable transformation of P.This indicates that if the adjoint state x db is given, the predicted value at the observation station can be evaluated just by calculating its spatial inner product with tsunami source x src .This is true regardless of bathymetry variation in either source or observation area, and nonpoint sources can be approximated via linear superposition of a suitable collection of point sources.The adjoint state x db can be precalculated and stored as a complete database if the location of the observation site is known beforehand, which is generally the case for coastal tide gauges or seafloor pressure gauges that record actual tsunami events.Note that the computational cost of the spatial inner product is much smaller than the time integration of the base equations.Furthermore, as shown above, the computation of this adjoint state database can be accomplished by performing simple variable transformations before and after the time integration of the original equations.Adjoint computation by variable transformation greatly reduces the cost of adjoint model development.Three ways of waveform evaluation at an observation site are shown in Figure 1.
The first way is based on the "forward" simulation with the base equations.Forward simulations were performed using a finite difference method with a staggered grid, and time evolution was calculated using the leapfrog method (cf., Baba et al., 2017;Goto et al., 1997).For the dispersion term, a Poisson equation for Φ is derived according to Shigihara and Fujima (2014), and its differentiated equations are solved by an iterative method.This part is the most computationally demanding in the simulation.To solve it efficiently, a parallel implementation of the conjugate gradient method of the Portable, Extensible Toolkit for Scientific computation (PETSc) library (Balay et al., 1997(Balay et al., , 2022) and a multigrid preconditioner, HYPRE (Falgout & Yang, 2002), is used.When including nonlinear terms in the simulation as a reference, the advection term was evaluated using first-order upwind differencing and the bottom friction term using central differencing, following Goto et al. (1997).The boundary with the shore was assumed to be perfectly reflective, and at the offshore boundary, an artificial "sponge" layer was placed outside the computational domain to absorb energy (Cruz, et al., 1993).Within the sponge layer, the dispersion term was gradually reduced toward the outer edge and was set to zero at the outer edge (Shigihara & Fujima, 2014).Radiation boundary conditions based on the long-wave approximation were then applied at the outer edge of the sponge layer (Goto et al., 1997).The deformation of the seafloor and changes of geoid are calculated by convolution of the tsunami load with the Green's function for the response due to a unit load (e.g., Pagiatakis, 1990;Ray, 1998;Vinogradova et al., 2015).Since this Green's function is axisymmetric, these can be efficiently computed by isotropic convolution using spherical harmonic transformation (e.g., Roddy & McEwen, 2021).A high-performance library Spherical Harmonic Transformation for numerical simulations (SHTns; Schaeffer, 2013) is used to perform this isotropic convolution.Gravitational acceleration is given by the Somigliana formula (Heiskanen et al., 1967).The detailed description of the numerical algorithm for the forward model is shown in Text S2 of Supporting Information S1.
The second way is based on "adjoint" simulation with the adjoint equations.This requires the adjoint model development.The adjoint model can be obtained by taking the transpose of the discretized forward model (cf., Kalnay, 2003;Appendix B).For the computational part of the deformation of the seafloor and changes of geoid, the self-adjoint property of isotropic spherical convolution (Kennedy et al., 2008) is utilized.The detail algorithm is shown in Text S2 of Supporting Information S1.
The third way is based on the "backward" simulation with the base equations and pre-and post-variable transformations.The fact that the backward and adjoint states can be converted by the transformation and inverse transformation according to Equation 20shows that the adjoint and backward simulations are consistent when Cartesian coordinates are used and the input and output are water levels, without considering density stratification, Boussinesq dispersion, and the Coriolis effect.The utilization of simple backpropagation without the variable transformation is justified only when all of these conditions are met.Equations 20-22 show that the adjoint and backward states are interchangeable through transformations that depend only on the density ratio of the seawater, the depth, and the latitude of the source and each grid point.This is why there is no phase difference between the adjoint and backward simulations but only a difference in amplitude in Figure 1.The third method brings the benefit of efficient waveform evaluation without the additional costs of adjoint model development.We call the latter two ways, which are numerically equivalent, the Adjoint Synthesis method.

SMC Method
A Markov Chain Monte Carlo (MCMC) method (e.g., Metropolis & Ulam, 1949) was used for parameter exploration.In the MCMC, the prior distribution p(m) of model parameters m and the likelihood p(d|m) are modeled and the posterior distribution p(m|d) can be asymptotically estimated by a series of random walks.
According to Bayes' theorem, the posterior distribution is explicitly expressed as follows: where k is a normalization constant.This constant is generally expressed as a multiple integral over parameter space and is extremely expensive to compute but does not affect the shape of the posterior distribution.MCMC can therefore estimate the posterior distribution from the results of a series of random walks without evaluation of the normalization constant.One of the most basic random walk methods is the Metropolis-Hastings algorithm (Hastings, 1970).In this method, a walker in the parameter space moves to a candidate location with the acceptance rate r defined as follows: where m cur and m can are the vectors representing the current and candidate locations of a walker in the parameter space, respectively.The locations of a series of random walks can be thought of as samples following a posterior distribution.In this way, we can estimate the posterior distribution.
The likelihood was defined as follows, assuming that the difference between the observed and modeled waveform variables follows a Gaussian distribution (e.g., Dettmer et al., 2016): where Σ is the covariance matrix, n is the number of waveform data d, and F is a function that converts model parameters to tsunami waveforms, which will be described in detail later.For simplicity, let Σ = σ 2 I, Here, σ 2 is the variance, and I is an identity matrix.Various distributions and their combinations can be used as a prior distribution.The case study shown later uses three distributions: uniform, normal, and exponential. Uniform: Normal: Exponential: where a and b are parameters that define the location and scale of the distribution, respectively.
For an earthquake, the posterior distribution of the fault parameters estimated using GNSS data is often multimodal (Ohno et al., 2021) and has a strong correlation between parameters (Munekane et al., 2016).Therefore, we use the SMC method (e.g., Buchholz et al., 2021;Nguyen et al., 2016), which enables efficient sampling from a multimodal distribution.
In the SMC method, a pseudo-inverse temperature β is introduced as follows.
This distribution coincides with the prior distribution when β = 0 and with posterior distribution when β = 1.A simple distribution, such as one of Equations 27-29, is typically chosen as the prior distribution, since they can be sampled using standard numerical libraries.The posterior distribution, on the other hand, can be a complex distribution with multiple peaks.The SMC method aims to approximate the posterior distribution with an ensemble of many samples.The scheme starts with the ensemble obtained at β = 0.As β is incrementally increased, the ensemble is adjusted to follow the distribution p β (m|d).Finally, an ensemble approximating the posterior distribution p(m|d) is obtained with β = 1.To be specific, the following "correction," "selection," and "mutation" are executed repeatedly.
Correction: Suppose that a set of samples {m i } and their weights {w i t 1 } that approximates p β t 1 (m|d) is known, where superscript i is the index of the sampled model, and subscript t 1 is the number of steps in the iteration of increasing β.The distribution is then expressed as follows: where N p is number of samples.The weights satisfy the following relationship.
Journal of Geophysical Research: Solid Earth 10.1029/2024JB028750 TAKAGAWA ET AL.
Due to the following relationship: the posterior distribution of β t can be expressed as follows: where By updating the weights according to the above equation, we can obtain the weights that approximate the posterior distribution of β t from the ensemble and weights that approximate the posterior distribution of β t-1 .
Selection: Continued updating of weights tends to result in one sample with large weights and many other samples with small weights, which causes degeneracy in approximation by ensembles.To avoid this, resampling is performed when the Effective Sample Size (ESS), represented by the following equation, falls below a threshold value.

ESS = ( ∑
Resampling is performed by sampling with replacement from {m i } with weights {w i t } .After the resampling, all new samples are given equal weights. Mutation: In the mutation step, the sample parameters are updated by the Metropolis-Hastings algorithm presented above.How the candidates are generated has a significant impact on the convergence of the algorithm.If the random walk step is taken too large, few will be accepted, and the samples will hardly move in the parameter space.On the other hand, if the steps are too small, the acceptance rate will increase, but the movement will be slower, and it takes more iterations to obtain an ensemble that well represents the target distribution.Here, random samples from a normal distribution are used as candidates (Nguyen et al., 2016).The mean μ and covariance Σ of the normal distribution are defined as follows: Journal of Geophysical Research: Solid Earth 10.1029/2024JB028750 TAKAGAWA ET AL.

Evaluation of Waveforms
The spatial distribution of slip on a rectangular fault is take into account in the function F of Equation 26, which transforms rectangular fault parameters into tsunami waveforms.This function is evaluated many times in the SMC algorithm.The slip distribution is inverted each time within this function.
Before describing the function F itself, how to evaluate tsunami waveforms for a basic uniform-slip rectangular fault model is presented.The procedure consists of the following four steps.(a) Calculate the horizontal and vertical displacement of the seafloor from fault parameters assuming a semi-infinite elastic medium based on the formulation by Okada (1992).(b) The vertical component of ocean displacement associated with the horizontal movement of the seafloor slope is added to the vertical displacement of the seafloor (Tanioka & Satake, 1996).(c) Calculate vertical displacement of the sea surface as a response of the seafloor displacement by using Kajiura (1963)'s formula.According to the formula, the water column acts as a depth filter with the factor of 1/ coshkh 0 in the wave number domain, where k is the wave number and h 0 is the water depth of the seafloor deformation area (e.g., Saito & Furumura, 2009;Ward, 2001b).(d) Calculate tsunami waveforms at observation sites by using adjoint synthesis.
To estimate the spatial distribution of slip, the rectangular fault is divided into N L segments in the strike direction and N W segments in the dip direction.The tsunami waveforms are calculated by giving a unit amount of slip to each of the N L × N W subfaults.These waveforms are used to estimate the amount of slip for each subfault by solving the following non-negative constrained least squares problem.
where x is a vector consisting of {x k }, x k is the slip on thekth subfault, the hat indicates that the value is an estimate, and d is a unified vector of observed waveforms at all observation points.σ is a standard deviation of the waveform modeling error.L is a Laplacian matrix such that Lx = 0 is equivalent to the following relationship: where i = {1,2,⋯,N L }, j = {1,2,⋯,N W }. L and W are the total length and width of the fault, respectively.α is a positive parameter that defines the spatial smoothness of the slip distribution.F is a matrix, whose k th column is a unified vector of modeled waveforms at all observation points from the k th subfault with unit slip.In the evaluation of the modeled waveforms, a time lag T gen between the earthquake occurrence time and the effective tsunami generation time and the modeling error T mod l of the tsunami propagation time at lth observation point are taken into account.Since the waveforms obtained by adjoint synthesis are discrete-time data, the time-shifted waveforms were obtained by linear interpolation.Some studies have considered the time variation of rupture for each sub-fault (e.g., Satake et al., 2013), but here, for simplicity, a common value is used for all faults as the time lag, and instantaneous fault rupture is assumed.(i.e., the effects of rise time and rupture propagation were not considered).According to Lotto et al. (2017), the effect of initial horizontal velocity on tsunami waveforms is negligible for tsunamis generated by typical subduction zone thrust fault ruptures.Initial horizontal velocity is assumed to be zero in this study.According to Tsai et al. (2013) and Watada et al. (2014), T mod l is expected to be smaller than 1% of the tsunami arrival time.In the method proposed here, these values of time lag of effective tsunami generation and modeling error of tsunami arrival time at each observation point are estimated in the SMC algorithm.
Equation 40 can be solved using the general algorithm of non-negative least squares, or NNLS (Lawson & Hanson, 1995), by rewriting it as follows: The algorithm proposed by Franc et al. (2005) is used in this study.The introduction of α, which controls the smoothness of the slip distribution, is intended to prevent the overfitting to the observed waveform data.To ensure that this α is well estimated in SMC, the idea of k-fold Cross Validation (CV) is used.In k-fold CV, the data are divided into k folds.One fold is selected and the remaining k 1 folds are used to estimate the slip distribution.The slip distribution is used to estimate data of the selected fold.This operation is repeated k times to obtain estimates for all folds.The estimates are used as the evaluation of F(m) in Equation 26.The entire SMC algorithm can be summarized as shown in Table 1.

2012 Haida Gwaii Earthquake and Tsunami
The 28 October 2012, Mw 7.8 earthquake off Haida Gwaii island, British Columbia, Canada (Figure 2), occurred as a result of shallow, oblique-thrust faulting near the plate boundary between the Pacific and North American plates (USGS, 2012).The earthquake source mechanism (e.g., Lay et al., 2013) and its tsunami impact (e.g., Fine et al., 2015;Leonard & Bednarski, 2014) have been extensively studied.The tsunami generated by this earthquake was observed at various locations over the Pacific Ocean, and waveforms were recorded by tide gauges, Deep-ocean Assessment and Reporting of Tsunamis (DART) buoy systems (González et al., 1998), and cabled bottom pressure recorders of Northeast Pacific Time-Series Undersea Networked Experiments (NEPTUNE)-Canada ocean observatory (Thomson et al., 2011).The fault parameters (Sheehan et al., 2015) and slip distribution (Gusman, Sheehan, et al., 2016) of the earthquake were estimated from these observed tsunami waveforms.

Tsunami Data and Numerical Models
The waveforms observed by 10 DART buoys were used (Figure 2).The original data with 15-s sampling were down-sampled to 30-s intervals and tidal components were removed by polynomial fitting (e.g., Wang & Satake, 2021).No filtering of any kind, such as bandpass filter, is applied.60-minute waveforms following the tsunami arrival were used to test the validity of the inverted source models.The ranges are indicated by the thick orange lines in Figure 2b.Vertical dashed lines indicate the post-event observation times at which tsunami source inversion analyses are considered below, so that we can consider how soon after the earthquake reliable results can be obtained.For the inversion, waveform data prior to each post-event observation time that fell within 15 min before and after the tsunami arrival (see Figure S1 in Supporting Information S1) was used.The exclusion of waveform data more than 15 min after the arrival was intended to avoid the effects of reflections and scattering from shallow waters.
For the tsunami simulation, bathymetry data with 2-arcminute and 30-arcsecond intervals were used.The data were generated by down-sampling the GEBCO_2019 grid to 15 arc-second intervals (GEBCO, 2019).The computational domain ranged between latitudes 64°north and 4°south, and the longitudes between 164°west and 60°west.The time step was set to 3 and 0.75 s for 2-arcminute and 30-arcsecond bathymetry data, respectively, to satisfy the Courant-Friedrichs-Lewy condition to avoid numerical instability (e.g., Goto et al., 1997).Five numerical models with different accuracies were used for fault parameter inversion and tsunami waveform estimation to examine the effect of the model accuracy on the estimation results.The effects considered by each model are summarized in Table 2. Model 1 is based on the linear long wave equation.Model 2 takes into account the effects of seawater density stratification, elastic loading and gravitational potential variation.These effects are collectively referred to below as secondary effects.The effect of Boussinesq dispersion is further incorporated in Model 3. Model 4 uses higher resolution bathymetry data with 30 arc-second grid intervals.Model 5 further takes into account nonlinear terms of advection and bottom friction.A value for Manning's roughness coefficient of 0.025 m 1 3 s, suitable for natural channels in good condition (Imamura et al., 2006), is used in Model 5.
The waveforms obtained at the 10 observation points using each of the five tsunami propagation models are shown in Figure 3.The rectangular fault with uniform slip proposed by Sheehan et al. (2015) was assumed as a tsunami source.The water depth of 2,000 m was chosen for the Kajiura (1963) filter as a representative value around the tsunami source area.Figure 3 shows the amount of time shift for each station.Observation points with smaller time shifts are located closer to the tsunami source.At DART46419, which is closest to the tsunami source, there is little difference between models 1 and 2. Comparison of models 2 and 3 shows that the onsets of the first tsunami waves are almost identical and the secondary effect is not remarkable.With respect to the first peak, Model 3 shows a slight time delay.This is due to Boussinesq dispersion, where the propagation velocity is smaller for the shorter period components.The impact of the secondary effects is clearly visible at distant observation points such as DART51407.Comparing models 1 and 2, the wave shapes are similar.The main difference is the tsunami arrival time.As Watada et al. (2014) showed, a small trough preceding the first peak is seen in DART 51407 due to the secondary effect.Tsai et al. (2013) estimated the reduction in tsunami velocity due to secondary effects compared to that predicted by the long wave approximation (i.e., our Model 1) for an average Pacific Ocean depth of 4,000 m.The effect of the compressibility of seawater is about 0.5%, while the effect of the elastic deformation of the earth depends on wavelength, and is 0.3% at a wavelength of 200 km.Watada et al. (2014) evaluated the effect of geopotential variation to be 0.06% for a period of 1,000 s, which corresponds to a wavelength of 200 km at a depth of 4,000 m.In total, the wave velocity reduction due to the secondary effects is about 0.86%.The Boussinesq dispersion term corresponding to this depth and wavelength reduces the wave speed by 0.27%, for a total reduction of approximately 1.1%.The difference between models 3 and 4 is due to the resolution of the bathymetry data used in the numerical simulation.The two waveforms are almost identical until the middle of the trough after the first peak, after which the difference becomes apparent.The effect of bathymetry data resolution is considered to be small for waves arriving directory from the tsunami source to the offshore observation points, while subsequent waves are more sensitive to the resolution because they have followed a more indirect path, reflecting off the coastlines or refracted by shallow bathymetry.The difference between models 4 and 5 is small, with only a slightly smaller amplitude in model 5 in the second hour of some waveforms.Again, this is thought to be mainly due to propagation through shallow water, where bottom friction more effectively attenuates the tsunami energy (Baba et al., 2017).

Adjoint State Database
Adjoint states for 9 hr were calculated using backward simulation with preand post-variable transformation for each observation point and for each model.Calculations were made using Models 1 to 4, but not Model 5, since its nonlinear terms violate the assumption of linearity on which adjoint synthesis is based.The computational domain, time step, and other conditions were the same as described in the previous section.Snapshots of adjoint states with 30-   s intervals were converted to grid data in a Cartesian coordinate system and stored as a database.The conversion was conducted by a transverse Mercator projection with center at 52.5°north and 132.25°west.The grid ranges ±150 km northward and eastward from this center, and the grid interval was 2 km.The reason for converting to a Cartesian coordinate system is to efficiently calculate seafloor displacement using Okada's (1992) formulae.This extent for the source area was set based on the W-phase moment tensor solution of USGS ( 2012) and the scaling relations of seismic moment, rupture area, and average slip proposed by Murotani et al. (2013).Reliable estimates of moment magnitude and source mechanism are assessed to be available within 25 min of the earthquake origin time for events larger than Mw 7.0 in real-time operations at USGS National Earthquake Information Center (Hayes et al., 2009).

Inversion Parameters
The geometry of a finite rectangular fault and slip distribution on it were estimated simultaneously from the observed tsunami waveforms by SMC.
The geometric parameters of the fault consist of northing and easting coordinates, depth D, strike, dip δ, rake, the aspect ratio defined as W/L, and the width ratio defined as Wsinδ/2D.The northing and easting coordinates and depth D are defined at the center of the fault.The width ratio ranges from 0 to 1.A width ratio of 1 implies a surface rupture, while a width ratio less than 1 implies an imbedded rupture (Savage & Hastie, 1966).In addition to these geometric parameters, the parameter α, which controls the smoothness of the slip distribution, observation error σ, effective tsunami generation time T gen , and the modeling errors of the tsunami propagation time T mod that were normalized by the arrival time of each observation point were estimated by SMC as samples from the posterior distribution defined by Equation 23.The number of samples was set to 256, and the intervals of pseudo-inverse temperature β t were determined adaptively based on a target ESS (Buchholz et al., 2021).Prior distributions were set as shown in Table 3.The zero-slip ratio in the table is the percentage of the number of sub-faults estimated to have zero slip.If the number of faults with zero slip exceeded 10%, the prior probability was set to zero, which we refer to below as the "zero-slip prior."Although SMC is an efficient method of sampling from a multimodal posterior distribution, it cannot completely eliminate the dependence on initial values.In order to minimize this dependence, sampling with SMC under the same conditions was conducted 20 times with different random seeds.Among the 20 ensembles, the ensemble containing the sample with the highest posterior probability was selected and used in subsequent analyses.
The slip distribution is solved by NNLS inside the SMC iteration.For CV, the waveforms used for inversion were divided into segments of 6 min each.Two folds of the odd-numbered segments and the even-numbered segments were used in the CV.After the sampling of SMC, the final slip distribution was inverted by NNLS again for the entire data set without folds.

SMC Sampling
The SMC sampling is expected to converge to the target posterior distribution as the pseudo-inverse temperature increases.Figure 4 shows the distribution of rectangular faults for each pseudo-inverse temperature in different colors.This figure shows the results for the Model 4 and the longest post-event observation time (henceforth we refer to this simply as "observation time") of 360 min.It can be seen that the extent of the fault converges between the Queen Charlotte Fault and Haida Gwaii Trench when the pseudo-inverse temperature is 1.This area is located  2015) source model of the 2012 Haida Gwaii earthquake, using tsunami computation Models 1 to 5. The tsunami arrival time at each observation point is shown in the figure.Observed waveforms are also shown as reference.The waveform of DART51407 is magnified four times.Model 1 is based on a basic linear long wave equation.Model 2 takes into account the secondary effects, that is, seawater density stratification, elastic loading, and geopotential variation.Model 3 takes into account these effects and Boussinesq dispersion.Model 4 is simulated using higher resolution bathymetry.Model 5 further considers nonlinear terms of advection and bottom friction.
Aspect ratio Uniform 0.5 0.2 Width ratio Uniform 0.5 0.5 a Probability density function.b Lay et al. (2013).on a topographic feature called the Queen Charlotte Terrace.This terrace is composed of a sedimentary wedge formed by the oblique subduction of the Pacific plate.This fault range indicates that the slip motion occurred between the Pacific Plate and the sedimentary wedge, and that the tsunamigenic slip was restricted to the seaward of the Queen Charlotte Fault.This result is consistent with other studies (e.g., Gusman, Sheehan, et al., 2016;Lay et al., 2013;Nykolaishen et al., 2015).

Journal of Geophysical
The posterior distribution of all parameters estimated by SMC are shown in Figure 5.This result is obtained using Model 4, with an observation time of 360 min.There are a variety of correlations between parameter pairs, including positive, negative and skewed correlations.This Bayesian approach allows us to quantify the uncertainty of each parameter considering not only the maximum likelihood value but also the correlation between the parameters.The first to third quartiles are shown in Figure 5a as a quantitative measure of uncertainty.Figure 5b illustrates the slip distribution of the fault for the sample that marked the best posterior probability.The posterior distribution of the amount of slip in each subfault and the modeling errors of tsunami propagation time for each observation point are shown in Figures 5c  and 5d, respectively.Many distributions shown in Figure 5a are asymmetric.
The flexibility to estimate such non-Gaussian posterior distributions is an advantage of Monte Carlo methods, including SMC.The estimated distributions can be used to quantitatively assess the uncertainty of the estimates.

Effects of Zero-Slip Prior and CV
Although SMC is an efficient method of sampling from a multimodal distribution, actual calculations with a limited number of samples cannot completely avoid being trapped by local peaks.This can cause the estimated value to change depending on the initial value.Here we show that the zero-slip prior and CV in SMC are effective in reducing initial value dependence.A comparison of the case with and without the zero-slip prior and CV in SMC is shown in Figure 6.The figure shows the results of SMC sampling runs with 20 different initial random values in each case.When the zero-slip prior was considered and CV was introduced in SMC (left panel), almost the same fault distribution was estimated regardless of initial values.When the zero-slip prior was not considered (center panel), or CV was not introduced (right panel), the estimated fault model differed significantly depending on the initial values.Therefore, the introduction of the zero-slip prior and CV proposed here is essential for obtaining stable estimates.

Effects of Propagation Model and Observation Time Length
The medians and interquartile range (IQR) of the posterior distribution of the initial sea surface height estimated by varying the propagation model and the observation time are shown in Figure 7. IQR is a measure of statistical dispersion, which is defined by the difference of third quartile Q 3 and first quartile Q 1 .By 60 min after the earthquake, the first wave of the tsunami was observed only at DART46419 (Figure 2), which is located south of the source.As the observation time increases, the number of useable observation waveforms increase, and the directional coverage also increases.By 90 min, the first peak was also observed at DART46404, located farther south.By 120 min, the first peaks were also observed at DART46407 more to south, and at DART46410 and 46409 to the west.After 120 min, the waveform data of the western DART buoys became available.By 180 min, the first waves were also observed at DART46411 more to the south and DART46403 more to the west.By 360 min, the first peaks were observed all 10 DART buoys, and a tsunami signal is available at DART51407.The tsunami source models of Sheehan et al. (2015) and Gusman, Sheehan, et al. (2016) are also presented for reference in Figure 7. Their models are characterized by an uplifted area along the Haida Gwaii Trench offshore of the epicenter.
In Model 1, the peak is estimated at the same location as the reference source models in the cases of 60-and 90min observation.If the observation time is longer than 120 min, the wave source is estimated much farther east.
Because Model 1 does not account for secondary effects, the modeled tsunami propagation speeds are larger than     2015) and Gusman, Sheehan, et al. (2016), respectively.The red star represents the relocated epicenter for the 2012 Haida Gwaii earthquake (Kao et al., 2015).
in Model 2 and the difference in arrival time grows larger for more distant stations (Figure 3).For longer observation times, the wave source is estimated at a greater distance from the observation points as a result of using model waveforms with a greater propagation speed bias.Yoshimoto et al. (2016) evaluated the impact of the secondary effects on tsunami source inversion for the 2010 Maule earthquake tsunami and showed a similar trend.
These results indicate that if the model is not sufficiently accurate, the accuracy of tsunami source estimation does not necessarily improve with increasing observational data.Model 2 estimates the peak at the same location as the reference source models up to 120 min because it accounts for secondary effects.However, in the 180-and 360min case, it still estimates the tsunami source at the wrong location.This suggests that the modeling error in Model 2 is smaller than in Model 1, but still large at distant observation sites.Model 3 further accounts for Boussinesq dispersion and estimated a sharper peak near the trench axis in the 120-min case.This is attributed to the improved reproduction of the propagation speed of waves with large wavenumbers generated by a sharp peak.In contrast to Model 2, for Model 3 there is no tendency for the accuracy of tsunami source location to decrease as observation time increases.The median has remained almost unchanged since 120 min, when the waveforms to the west became available in addition to the ones of southern stations.The increase in observation time contributes rather to the reduction of uncertainty.Model 4, using higher-resolution bathymetry, estimated a tsunami source with a peak height close to that of the final tsunami source model, even at the shorter 60-minite observation time.
Comparing the waveforms of DART46419 for Models 3 and 4 within 60 min in Figure 3, there is only a very slight difference in the phase of the second peak.These small differences may have contributed to the improved accuracy of tsunami source estimation with waveforms of short (60 min) duration.
Comparisons of the observed and estimated waveforms are shown in Figure 8.The estimated waveforms were computed based on the parameters that showed the maximum posterior probabilities in the ensemble SMC solution.The durations of the waveforms displayed there are the same as the orange lines in Figure 2.Each column shows the results of a different propagation model.Waveforms based on tsunami source models proposed by Sheehan et al. (2015) and Gusman, Sheehan, et al. (2016) are also shown for reference.Their waveforms were calculated using Model 4. Color indicates the post-event observation time used for the inversion with SMC, with the interval of waveform data used for each observation time drawn with lines of the same color as the corresponding waveforms.The coefficient of determination R 2 is used in this study as a measure of waveform estimation accuracy.R 2 is defined as follows: where y and ŷ denote observed and modeled data, and y is the mean of the observed data.The value of R 2 is calculated using waveforms within the interval of 60 min shown in Figures 2 and 8, and the arithmetic mean of all observation points is also shown as the overall mean.The effective tsunami generation time and propagation time modeling errors are taken into account to estimate the waveform and calculate R 2 .To account for modeling errors, the R 2 of the published source models, were taken to be the maximum value when the time shift was varied from 3 to 3 min.
For Model 1, the R 2 is low, and accuracy does not improve as observation time increases.In particular, the estimation accuracy at the three distant observation points of DART 46413, 46402, and 51407 is poor.Although the estimated waveforms in the inverted range at the DART buoys nearer the source have good fits, the discrepancy in the later parts of these waveforms are large.This may be due to the effect of incorrectly estimated tsunami sources (Figure 7).Model 2 takes into account the secondary effects, which clearly improves the accuracy of the three distant locations.However, accuracy sharply degrades when the observation time exceeds 120 min.This can be attributed to the inadequate accuracy of the propagation model and the incorrect location of the estimated tsunami sources (Figure 7).Model 3 takes into account the Boussinesq dispersion, which is expected to improve the reproducibility of wave shapes, especially at distant observation points.In fact, the accuracy of the most distant point of DART51407 was significantly improved compared to that of Model 1 and 2. Furthermore, it differs from Models 1 and 2 in that the accuracy increases as the observation time increases.Model 4 uses high-resolution bathymetry to reproduce the waveforms even more accurately than Model 3. The difference is especially noticeable at DART 51407.The difference between the sources estimated in Models 3 and Where no inverted range of that color appears, it means the tsunami has not arrived at the DART buoy by that observation time, so there is no data to invert.Note that the latter half is not used for inversion, even if the observation time is long.The data for DART51407 is enlarged by a factor of 2 in height.The waveforms calculated by Model 4 from the source models of Sheehan et al. (2015) and Gusman, Sheehan, et al. (2016) are also shown as Models S and G for reference.
The bar charts on the right show the R 2 scores of these models and observation times.
4 is not large, being mainly evident at DART 51407, suggesting that the benefit of the high-resolution bathymetry was important for tsunami energy passing through the complex bathymetry around the Hawaiian Islands.On average of all DART buoys, R 2 , a measure of waveform prediction accuracy, is greater than 70% of the final fault model in the case of 60/90-min post-event observation.This suggests that highly accurate predictions are possible even from a very limited set of waveform data.The largest improvement in accuracy occurs when the post-event observation time increases from 90 to 120 min.Up to 90 min, the tsunami waveforms were captured only at stations in the south direction, but at 120 min, data from two stations in the west were added.The first wave observed at DART 46410 has a broad peak, while DART46409 records a characteristic waveform with a double peak.The improved directional coverage and observation of characteristic waveforms are considered to have resulted in a relatively large improvement in accuracy at 120 min.The distribution of tsunami sources was elevated along the trench in the estimation up to 90 min, but changed to a slightly oblique distribution to the trench after these waveforms were observed.
The effects of propagation model and post-event observation time on the estimation of tsunami sources and waveforms can be summarized as follows.Secondary effects were important for source location estimation, even for waveforms observed within 90 min observation time.Adding Boussinesq dispersion to these effects further improved source estimation when inverting observations from distant locations.Using in addition high-resolution bathymetry data further improved the accuracy of tsunami source estimation, even for observation times as short as 60 min.It also improved the accuracy of waveform estimation, especially for waveforms passing through shallow and complex bathymetry.The combination of highly accurate propagation models and exhaustive parameter exploration with the help of adjoint synthesis and SMC provides three major advantages: (a) high accuracy of tsunami waveform estimation, especially in the far field, (b) accuracy of waveform estimation increases monotonically as observation data increases, and (c) convergence to a final accurate solution in a relatively short observation time.In particular, the ability to predict distant tsunami waveforms with high accuracy from a limited amount of observation data near the tsunami source is a potentially important advance in terms of early warning for distant areas.

Uncertainty of Waveform Estimation
The SMC provides an ensemble of estimates of the tsunami source.An ensemble of estimated waveforms can also be obtained by applying the adjoint synthesis method to each tsunami source in the ensemble.The range of uncertainties in the tsunami waveform estimated from the ensemble in Model 4 is shown in Figure 9.The figure shows the interquartile range (IQR) and a measure of confidence range, which is defined the range from Q 1 1.5 IQR to Q 3 + 1.5 IQR (Tukey, 1977).In the case of the 60-min observation, only the tsunami waveform at one southern station DART46419 is available.The median of the estimated waveform did not reproduce the doublepeaked first wave feature that is prominent in western DARTs 46402, 46403, and 46409.However, the uncertainties were estimated to be relatively large, and most of the observed waveforms were in the confidence interval, except for the later peaks seen around 40 min in southern DARTs 46419, 46404, and 46407.In the 120min case, waveforms at two western DARTs are available.The improved directional coverage of observation sites and the characteristic double-peaked waveform data further improved the reproducibility of the waveforms at the farther western DARTs 46403, 46402, and 46413 and of the ebb waveforms around 35 min at southern DARTs 46419, 46404, and 46407, which were not used in the inversion.The later peaks at DARTs 46419, 46404, and 46407 were still out of the uncertainty range, but sharper peaks were estimated compared to the 60min case, showing an improvement in accuracy.In the 360-min case, changes from the 120-min case were relatively minor, but further improvement in waveform reproducibility and reduction in confidence intervals are observed.

The Final Fault Model
The results of our final finite fault model estimated by SMC sampling with Model 4 and data available at 360-min observation time are summarized in Figure 10 and its parameters are shown in Table S1.The areas of large slip are distributed in a linear fashion on the fault plane (Figure 10a), and the direction is slightly counterclockwise from the fault strike direction.The comparison of the displacement estimated by Okada's formula from the final model with the displacement observed by land GPS stations (Nykolaishen et al., 2015) shows good agreement, even though the inversion uses only tsunami waveform data (Figures 10b and 10c).Observed waveforms and those estimated using Model 4, which can be evaluated quickly using adjoint synthesis, and Model 5, which includes nonlinear terms of advection and bottom friction and requires forward simulation, are shown in Figure 10d.Models 4 and 5 both fit the observations well, not only the interval used for inversion, but also for the later peaks as indicated by the black arrows.There are little differences between Models 4 and 5 in the first waves, but about an hour after the tsunami arrivals, the amplitudes of Model 5 become slightly smaller, and the results are closer to the amplitude of the observed waveforms.The later waves contain a large component of waves that have passed through shallow waters.The energy of these components is attenuated by bottom friction, and their amplitude is reduced (e.g., Baba et al., 2017;Saito et al., 2014).Although Model 4 tends to overestimate the amplitude slightly in the later waves due to the lack of bottom friction, it presents waveform predictions with good accuracy over a long period of time.The waveforms of Model 4 can be computed in a short time by adjoint synthesis and are expected to provide reliable information for making decisions on the issuance or cancellation of tsunami warnings, especially for distant areas.It is necessary to consider additional detailed bathymetric and topographic effects including inundation when assessing tsunami impact in shallow water and on land.To evaluate these effects, the estimated tsunami source or offshore waveforms could be used as input to nonlinear tsunami propagation models (e.g., Oishi et al., 2015;Titov et al., 2016).
The scope of application of the adjoint synthesis method is not limited to the fault parameter estimation shown here.Since tsunami waveforms generated by arbitrary nonpoint sources can be evaluated by superposition as long as the linear approximation is valid, it is possible to extend the analysis to the case of multiple faults and curved faults, and even to apply the method to tsunamis generated by moving sources, such as submarine landslide tsunamis (e.g., Ward, 2001a), meteotsunamis (Monserrat et al., 2006), and tsunamis generated by atmospheric Lamb waves (Kubota, et al., 2022).In these cases, the advantages of adjoint synthesis can be utilized because the situation is the same as in the above examples: the location of the observation point is known before the event and the number of observation points is relatively small, compared to the extremely high degree of freedom of the tsunami source.It can also be applied to Green's Function-based Tsunami Data Assimilation (GFTDA) proposed by Wang et al. (2017).GFTDA is a method used to speed up the data assimilation of tsunami (Maeda et al., 2015) using a pre-computed database, but the extensive computation time required to compute the database can be prohibitive.The database of densely gridded Green's Functions can be computed using adjoint synthesis with much less cost in terms of time, hardware and electricity.The adjoint model and variable transformation methods proposed in this study enable adjoint operations that take into account the secondary effects, Boussinesq dispersion, and Coriolis effect in a spherical coordinate system, that are essential for high accuracy in transoceanic tsunami simulation.
The data range used for inversion in this study was 15 min before and after the tsunami arrival, but later phase data could be used for faster and more accurate tsunami source and waveform estimation.However, the use of later waveforms, which have relatively larger modeling errors, may adversely affect the accuracy of the inversion.What intervals should be used for inversion is not examined here.The applicability of the proposed method to tide gauge data, which is sensitive to shallow water bathymetry and bottom friction, has not been investigated.Although only data from the distant DART buoys, which are less affected by seismic waves, were used, it may be possible to predict tsunamis earlier if data from ocean bottom pressure gauges near the tsunami source can be utilized.On the other hand, noise from seismic motions and crustal deformation may cause problems.How to handle and integrate multiple types of observation data with different nature of errors, such as modeling errors, boundary condition errors, and observation errors, is an important issue left for future research.

Conclusions
We developed a method for the synthesis of tsunami waveforms at observation sites generated by arbitrary nonpoint tsunami sources, which utilizes a precomputed adjoint state database.This database does not require any assumptions about the tsunami source, unlike other forward model databases which use a priori assumptions about earthquake generation based on specific fault geometries.The database itself can be efficiently computed by using adjoint operators, because such computations scale by the number of observation sites instead of by the typically much larger set of points used to represent potential sources.The adjoint operation on accurate governing equations for distant tsunamis is shown to be easily implemented by a typical forward computation with pre-and post-variable transformations.
The adjoint synthesis method was applied to tsunami waveform inversion for simultaneous estimation of fault geometry and slip distribution.The fault geometry was estimated by SMC.The slip distribution was estimated by non-negative least squares with regularization of spatial smoothness in the iterative steps of the SMC.CV was introduced into the waveform evaluation step in the SMC to suppress overfitting.It was found that the convergence to the global optimal solution is enhanced by placing an upper bound on the zero-slip ratio as part of the prior probability.
The proposed method was applied to the 2012 Haida Gwaii earthquake tsunami to estimate the tsunami source and waveforms.It was examined how secondary effects, that is, density stratification, deformation of the seafloor and geoid due to tsunami loading, and Boussinesq dispersion, fineness of bathymetric data, and nonlinear terms of advection and bottom friction affect estimation of tsunami source and waveforms.By taking into account the secondary effects, the reproduction accuracy of the tsunami arrival time at distant stations was increased, and the estimation accuracy of the tsunami source location was improved.When Boussinesq dispersion was considered, it was estimated that slip spread in a linear fashion oblique to the strike.
In order to understand the potential of using this approach in a tsunami warning system, the effect of the postevent observation time on the estimation of the tsunami source and waveforms was investigated, and it was found using less accurate propagation models causes the accuracy of the tsunami source estimation to decreases as the observation time increases, contrary to what might be expected.In contrast, a highly accurate model that took into account secondary effects, Boussinesq dispersion, and detailed bathymetry resulted in an estimate source accuracy that increases with increasing observation time.This model was able to estimate the approximate distribution of tsunami sources based on 60 min of waveform data with half-waves of the first wave observed at only one station to the south.In the 120-min observation case, tsunamis were observed at two western stations, improving the directional coverage, and a tsunami source with an uplifted area extending obliquely from the trench axis, improving the accuracy of waveform estimation.When the observation time was extended beyond 120 min, the median value of the estimated tsunami source did not change significantly, but the uncertainty gradually decreased.
The Final Fault Model, which was estimated using a highly accurate propagation model, Model 4, and the longest waveform data of 360 min reproduced the displacement of the GPS stations on land well.Tsunami simulation was also performed with the consideration of nonlinear terms of advection and bottom friction.The nonlinear tsunami simulation shows almost the same waveform as the linear simulation for about 30-60 min after the first tsunami arrival, but the amplitude becomes slightly smaller in late-arriving energy, indicating a better correspondence with the observed waveform amplitude changes.Although the influence of the nonlinear terms cannot be considered in the adjoint synthesis method, its influence is limited to the waveform in the latter parts of the waveforms and even then, the amplitude difference is quite minor.The inversion proposed here, which combines the adjoint synthesis method and SMC, can efficiently estimate the tsunami source and waveforms and their uncertainties quantitatively.The quantitative prediction would help in the proper updating and canceling of tsunami warnings.

Figure 1 .
Figure1.Three ways of waveform evaluation at an observation point.The usual forward model simulates tsunami propagation from the tsunami source and the water level is evaluated at the observation point.In the adjoint synthesis method, on the other hand, the waveform is evaluated by placing a point tsunami source at the observation point, simulating propagation using an adjoint model, and taking the spatial inner product with the tsunami source used in the forward model.The adjoint model simulation can be replaced by the backward simulation with the usual forward model using pre-and post-computation variable transformations.The three waveforms are shown with intervals in the center panel, each obtained by these different ways of waveform evaluation.These waveforms show excellent match even with a nonpoint tsunami source over a variable bathymetry.

Figure 2 .
Figure 2. (a) Focal mechanism of the 2012 Haida Gwaii earthquake and observation points of the DART buoy system.White lines indicate the contour of tsunami arrival time.The black lines represent 5-mm contours of the maximum tsunami height.(b) Tsunami waveforms observed by DART buoys.The orange parts indicate the interval used to verify statistical accuracy.Dashed vertical lines indicate the post-event observation times (60, 90, 120, 180, and 360 min) at which source inversion analysis is undertaken in Section 4.

Figure 3 .
Figure 3.Comparison of waveforms at 10 observation points simulated for theSheehan et al. (2015) source model of the 2012 Haida Gwaii earthquake, using tsunami computation Models 1 to 5. The tsunami arrival time at each observation point is shown in the figure.Observed waveforms are also shown as reference.The waveform of DART51407 is magnified four times.Model 1 is based on a basic linear long wave equation.Model 2 takes into account the secondary effects, that is, seawater density stratification, elastic loading, and geopotential variation.Model 3 takes into account these effects and Boussinesq dispersion.Model 4 is simulated using higher resolution bathymetry.Model 5 further considers nonlinear terms of advection and bottom friction.

Figure 4 .
Figure 4. History of convergence of rectangular faults by sequential Monte Carlo sampling.Fault rectangles are drawn in different colors for each pseudo-inverse temperature.

Figure 5 .
Figure 5. (a) Corner plot of the fault geometric parameters obtained using the sequential Monte Carlo method.The contour lines represent first, second and third quartiles.Q 1 , Q 2 , and Q 3 , denote first, second and third quartiles, respectively.(b) Slip distribution of the sample with the best posterior probability.(c) Posterior distribution of slip amount of each subfault.(d) Posterior distribution of modeling error of tsunami arrival time, which is normalized by tsunami arrival time.

Figure 6 .
Figure 6.Fault models estimated from 20 different initial values.(left) proposed method.(center) Cases without zero-slip prior.(right) Cases without Cross Validation.

Figure 7 .
Figure 7. Median and interquartile range (IQR) of the posterior ensemble of initial water level distributions estimated by varying the propagation model and observation time.Model S and G indicates the published source models proposed by Sheehan et al. (2015) andGusman, Sheehan, et al. (2016), respectively.The red star represents the relocated epicenter for the 2012 Haida Gwaii earthquake(Kao et al., 2015).

Figure 8 .
Figure 8. Observed tsunami waveforms and estimated waveforms by the four propagation models with varying post-event observation times.The arrival time at each observation point is shown below the DART buoy number.The relative time section between 15 and 45 min from arrival shown here is the same as the range shown by the orange lines in Figure 2b.For better clarity, these waveforms are plotted against the time relative to the tsunami arrival.In Figure S1 of Supporting Information S1, they are also plotted against absolute elapsed time to see the overall time.Estimated waveforms are colored according to the observation time indicated in the legend (60-360 min).Where no inverted range of that color appears, it means the tsunami has not arrived at the DART buoy by that observation time, so there is no data to invert.Note that the latter half is not used for inversion, even if the observation time is long.The data for DART51407 is enlarged by a factor of 2 in height.The waveforms calculated by Model 4 from the source models ofSheehan et al. (2015) andGusman, Sheehan, et al. (2016)  are also shown as Models S and G for reference.The bar charts on the right show the R 2 scores of these models and observation times.

Figure 9 .
Figure 9. Observed waveforms and the ranges of uncertainty in waveforms estimated at different post-event observation times for Model 4. Uncertainty is indicated by a confidence range based on IQR.The time interval used for inversion is indicated by green lines.For clarity, these waveforms are plotted against the time relative to the tsunami arrival.The arrival time at each observation point is shown below the DART buoy number.These waveforms are also plotted against absolute elapsed time to see the overall time in Figure S2 of Supporting Information S1.

Figure 10 .
Figure 10.(a) Slip distribution of the final model.The gray arrow indicates the direction of slip on the sub-fault where maximum slip was estimated.(b and c) Horizontal and vertical displacements at GPS stations.Observation data is from Nykolaishen et al. (2015).(d) Observed tsunami waveforms and estimated waveforms by Models 4 and 5.

Table 1
Algorithm of Simultaneous Estimation of Fault Geometric Parameters and Slip Distribution can } from normal distribution defined by {μ} and {Σ} Compute horizontal and vertical displacement at sea bottom generated by a unit slip of each subfault defined by {m can } Compute total vertical displacement including slope effect Compute vertical displacement of sea surface by Kajiura filter Compute waveforms at observation points by adjoint synthesis Compute {F(m can )} by solving NNLS with k-fold CV Compute acceptance rate r Update {m} by the Metropolis-Hastings algorithm End for End for Compute horizontal and vertical displacement at sea bottom generated by a unit slip of each subfault defined by {m} Compute total vertical displacement including slope effect Compute vertical displacement of sea surface by Kajiura filter Compute waveforms at observation points by adjoint synthesis Compute {F(m)} and slip amounts by solving NNLS with all waveform data Note.{}indicates that it is a set of samples.TAKAGAWA ET AL.

Table 2
Five Models and What Is Considered a Seawater density stratification, seafloor deformation due to tsunami loading, and gravitational potential change.b Finer grid resolution of 30 arcseconds otherwise grid resolution is 2 arcminutes.c Advection and bottom friction.

Table 3
Prior Distribution of Parameters Used in Sequential Monte Carlo