Global Dissipation Models for Simulating Tsunamis at Far-Field Coasts up to 60 hours Post-Earthquake: Multi-Site Tests in Australia

At far-field coasts the largest tsunami waves may occur many hours post-arrival, and hazardous waves may persist for more than 1 day. Such tsunamis are often simulated by nesting high-resolution nonlinear shallow water models (covering sites of interest) within low-resolution reduced-physics global-scale models (to efficiently simulate propagation). These global models often ignore friction and are mathematically energy conservative, so in theory the modeled tsunami will persist indefinitely. In contrast, real tsunamis exhibit slow dissipation at the global-scale with an energy e-folding time of approximately 1 day. How strongly do these global-scale approximations affect nearshore tsunamis simulated at far-field coasts? To investigate this we compare modeled and observed tsunamis at sixteen nearshore tide-gauges in Australia, generated by the following earthquakes: M w 9.5 Chile 1960; M w 9.2 Sumatra 2004; M w 8.8 Chile 2010; M w 9.1 Tohoku 2011; and M w 8.3 Chile 2015. Each tsunami is represented using multiple published source models, to prevent bias in any single source from dominating the results. Each tsunami is simulated for 60 h with a nested global-to-local model. On nearshore grids we solve the nonlinear shallow water equations with Manning-friction, while on the global grid we test three reduced-physics propagation models which combine the linear shallow water equations with alternative treatments of friction: 1) frictionless; 2) nonlinear Manning-friction; and 3) constant linear-friction. Compared with data, the frictionless global model well simulates nearshore tsunami maxima for ≃ 8 h after tsunami arrival, and Manning-friction gives similar predictions in this period. Constant linear-friction underestimates the size of early arriving waves. As the simulation duration is increased from 36 to 60 h, the frictionless model increasingly overestimates observed wave heights, whereas models with global-scale friction work relatively well. The constant linear-friction model can be improved using delayed-linear-friction, where propagation is simulated with an initial frictionless period (12 h herein). This prevents systematic underestimation of early wave heights. While nonlinear Manning-friction offers comparably good performance, a practical advantage of the linear-friction models herein is that solutions can be computed, to high accuracy, via a simple transformation of frictionless solutions. This offers a pragmatic approach to improving unit-source based global tsunami simulations at late times.


INTRODUCTION
Large earthquake-tsunamis can be hazardous even at far-field coasts. For example the 1946 Aleutian earthquake-tsunami caused 159 deaths in Hawaii; the 1960 Chile earthquaketsunami caused 142 deaths in Japan; the 2004 Sumatra earthquake-tsunami led to 300 deaths in Somalia (Fritz and Borrero, 2006;Okal, 2011). Smaller yet more-common tsunamis are also of interest for risk management because they can generate hazardous currents and minor inundation, with potential to harm people, damage assets, and disrupt economic activity at ports (e.g., Beccari, 2009;Borrero et al., 2015b). To mitigate the risk, early warning systems and hazard assessments are employed to guide emergency management and planning (Kânoglu et al., 2015;Grezio et al., 2017;Davies and Griffin, 2018). They are underpinned by numerical models of tsunami propagation and inundation which exploit various approximations for computational tractability (e.g., Gica et al., 2008;Greenslade et al., 2011;Lorito et al., 2015;Setiyono et al., 2017;Volpe et al., 2019;Davies and Griffin, 2020). It is important to understand how well such models simulate key quantities of interest for applications (e.g., maximum wave heights and current speeds, their timing, and the duration of dangerous waves). Herein we focus on the tsunami wave size at far-field coasts in the period from 0-to-60 h post-earthquake. This is motivated by the fact that hazardous waves can persist for several days (Tang et al., 2012). Some historical tsunami observations in Australia and New Zealand exhibited wave maxima arriving 1-to-2 days post earthquake, long after the initial tsunami arrival (Pattiaratchi and Wijeratne, 2009;Borrero et al., 2015a).
For tsunami modeling applications at far-field sites, it is often advantageous to simulate global-scale propagation with linear models, noting nonlinearity should be small because ocean depths are generally much larger than tsunami amplitudes (Shuto, 1991). Linear models are typically faster to solve numerically (e.g., Liu et al., 2008;Baba et al., 2014), and if many scenarios are required then even greater speedups follow by constructing solutions S(x, t) as linear combinations "unitsource" solutions U i (x, t).
S (x, t) i ∈ set of unit−sources in database Here x, t denote space and time, while the s i are constant coefficients defining the solution S (e.g., Gica et al., 2008;Miranda et al., 2014). For linear models S will be an exact solution, because the definition of linearity implies linear combinations of solutions are also solutions, and the calculation is very fast once a unit-source database is constructed (containing solutions U i ). For this reason the unitsource approach is very popular for large-scale probabilistic hazard assessment (e.g., Burbidge et al., 2008;Li et al., 2016;Molinari et al., 2016;Davies et al., 2017;Davies and Griffin, 2020;Zhang and Niu, 2020). Unit-sources also greatly simplify tsunami source-inversion algorithms, including for early-warning applications, because techniques from linear-regression can be combined with tsunami observations to solve for s i (e.g., Tang et al., 2012;Fujii and Satake, 2013;Percival et al., 2014;Romano et al., 2016). Unit-source solutions U i are sometimes computed with nonlinear hydrodynamic models (e.g., Yue et al., 2015;Molinari et al., 2016;Zhang and Niu, 2020) in which case Eq. 1 does not produce an exact solution of the original model, but irrespective solutions derived from linear combinations of unitsources are linear (by construction). In practice this approach works well if nonlinearity is small, which can be tested on a caseby-case basis (e.g., Yue et al., 2015;Molinari et al., 2016;Zhang and Niu, 2020). Linear models cannot simulate inundation and become unreliable if the wave amplitude is a significant fraction of the water depth, but are widely used to force nonlinear coastal inundation models using one-way or two-way nesting (e.g., Tang et al., 2009;Baba et al., 2014;Borrero et al., 2015a).
Tsunami propagation models at ocean-basin scales also often neglect friction, which has little effect on deep ocean tsunami simulations for a few hours post-arrival (Tang et al., 2012;Fujii and Satake, 2013;Allgeyer and Cummins, 2014;Baba et al., 2017;Heidarzadeh et al., 2018;Davies, 2019). However the frictionless shallow water equations are energy conservative, assuming smooth solutions, which mathematically implies the tsunami persists forever (Arakawa and Hsu, 1990;Fjordholm et al., 2011;Tang et al., 2012). In contrast real global-scale tsunamis eventually dissipate. Observations at coastal and deep-ocean gauges in the Pacific and Indian Oceans suggest an exponential time-decay of available potential energy, with an e-folding timescale about 1 ± 0.5 days following the initial diffusion of the tsunami energy throughout the ocean (Miller et al., 1962;Munk, 1963;van Dorn, 1984, van Dorn, 1987Mofjeld et al., 2000;Rabinovich et al., 2011;Nyland and Huang, 2013;Rabinovich et al., 2013). The wave-amplitude e-folding timescale is twice as long, because available potential energy is proportional to the squared wave-amplitude. Empirically, the e-folding timescale depends on the tsunami spectral properties as well as the site and its distance from the source (Rabinovich et al., 2011;Rabinovich et al., 2013). Comparatively short energy e-folding timescales have been reported in small coastal seas (x3 − 13 hours; van Dorn, 1987;Oh and Rabinovich, 1994) which can be explained if tsunami dissipation mainly occurs in shallow shelf regions (van Dorn, 1987). Munk (1963) suggested that tsunami dissipation might additionally be affected by energy transfer to internal ocean waves. This was subsequently confirmed to be important for simulating tidal dissipation (Munk, 1997;Egbert and Ray, 2000;Llewellyn Smith and Young, 2002) and is typically represented in tidal models using spatially varying linear friction, in addition to standard quadratic bottom-friction (Buijsman et al., 2015;Kleermaeker et al., 2019), although the significance of this process for tsunami dissipation remains unclear.
Quadratic bottom friction (e.g., Manning or Chezy) is most common in tsunami models, but the solutions cannot be represented with unit-sources (Eq. 1) over timescales for which the nonlinear dissipation is important. To represent global-scale tsunami dissipation within a linear hydrodynamic model, Fine et al. (2013) and Kulikov et al. (2014) combined the linear shallow water equations (LSWE) with a constant linearfriction term. Although not justified from hydrodynamic theory, this model implies an exponential time-decay of the tsunami's energy which is consistent with observational studies. The linearfriction coefficient may thus be estimated using observed tsunami decay timescales (Fine et al., 2013;Kulikov et al., 2014). The model was implemented using numerical methods that have good energy conservation to prevent numerical dissipation from dominating the results at late times, which was reportedly easier to achieve by modifying a linear model (Fine et al., 2013;Kulikov et al., 2014). Numerical dissipation is often significant for models solving the nonlinear shallow water equations (NSWE) and can even exceed the physical dissipation of global scale tsunamis at computationally practical resolutions (Popinet, 2011;Tang et al., 2012;Tolkova, 2014). This has been used to explain persistent under-estimation of late time tsunami wave heights in some nested nearshore models, even without any friction in the global propagation model (Tang et al., 2012;Tolkova, 2014). However, in the absence of significant numerical dissipation, friction is necessary to represent the late-time energy decay observed in global-scale tsunamis.
Modeled coastal tsunami wave heights are likely to be qualitatively affected at sufficiently late times by the treatment of global-scale dissipation. This has the potential to affect tsunami hazard assessments (e.g., runup maxima) and warnings (e.g., warning cancellation). Our study seeks to better understand the practical significance of this issue, focusing on the empirical performance of alternative models rather than the physics of tsunami dissipation. To this end we test several global-to-local scale nested-grid tsunami models by comparison with tsunami observations at tide-gauges in Australia, for the period 0-60 h post-earthquake. The alternative models differ only in their treatment of friction on the global-scale grid (where propagation is simulated with some variant of the LSWE); in all cases the tsunami at nearshore sites of interest is simulated on nested grids which solve the NSWE with Manning-friction. The tsunami initial conditions are derived from published earthquake-source inversions. We focus on relatively simple tsunami models which are computationally cheap and very practical in applications, while acknowledging that tsunami propagation may be simulated more accurately using computationally intensive approaches that include higher order physics such as dispersion, loading, seawater density stratification, and self-gravitation (Allgeyer and Cummins, 2014;Baba et al., 2017). Model improvements may also be sought through use of higher quality elevation data outside the primary area of interest, combined with higher model resolution, to better represent remote reflections (Kowalik et al., 2008;Geist, 2009). The practical benefit of very high-accuracy propagation modeling may be limited for hazard assessments, because plausible source variations can have an even greater effect on the simulated tsunami (Davies, 2019), although further study of this issue is warranted.
The paper is structured as follows. Section 2.1 presents the earthquake-tsunami sources (initial water-surface perturbations) used to model each historic event. Section 2.2 presents the tidegauge data used to test each model. Section 2.3 reviews the tsunami model setup. Section 2.4 details the alternative reducedphysics hydrodynamic models that are tested. Section 3.1 illustrates each model's performance using examples from three historic tsunamis. Sections 3.2 and 3.3 compare the modeled and observed tsunami maxima for all sites and events. Section 3.4 considers how the model errors are affected by the tsunami source.

Earthquake Source Models for Historical Tsunamis
Five historic tsunamigenic earthquake events were analyzed in this study: M w 9.5 Chile 1960, M w 9.2 Sumatra 2004, M w 8.8 Chile 2010, M w 9.1 Tohoku 2011, and M w 8.3 Chile 2015 ( Figure 1). These events were chosen because they generated relatively large tsunamis and resulted in good quality tide-gauge observations in Australia. The earthquake sources were represented using twelve published finite-fault inversions ( Figure 2) which represent an ad-hoc sample from the literature. The set of finite-fault inversions was not modified on the basis of our tsunami model performance because that could inadvertently compensate for any model biases which are of primary interest herein. Instead multiple inversions were used for each historic tsunami to prevent bias in a single source model from dominating the results.
Most tsunami initial conditions in Figure 2 were derived by computing the vertical co-seismic deformation from the published fault-geometries and slip vectors, assuming instantaneous rupture in a homogeneous elastic half-space (Okada, 1985;Meade, 2007). The effect of horizontal deformation over a sloping seabed (Tanioka and Satake, 1996) was not included as this tends to add shorter waves to the source which may not be well simulated with our long-wave hydrodynamic model (Saito et al., 2014). However source R14 was taken directly from the finite-element model of Romano et al. (2014); this includes horizontal components so is rougher than our other sources ( Figure 2). Source H19 was derived from watersurface unit-sources following Ho et al. (2019). A Kajiura filter (Kajiura, 1963;Glimsdal et al., 2013) was applied to all models except the water-surface inversion H19.
The studies on which our source models are based generally inverted the source using deep sea and/or coastal tsunami data, sometimes combined with geodetic data (GPS, InSAR, leveling) and occasionally teleseismic data. None made use of the Australian tide-gauge data studied herein. All inverted the source with linear combinations of earthquake and/or tsunami Green's functions, in some cases including regularization constraints (L11, Y14, R16, W17, Y18). Most simulated the co-seismic displacement assuming a homogeneous elastic halfspace, although R14 used a 3D finite-element model to account for material heterogeneity. Tsunami Green's functions were mostly modeled with shallow water schemes, although R14, R16 and Y18 employed a non-hydrostatic model. Some papers reported both a best-fitting source model and an average source model; in such cases we always used the former. Full details are available in the repository (see Data Availability Statement).

Test Sites and Data
The model is compared with data at 16 sites in Australia that have good bathymetry and relatively good quality tsunami observations at tide-gauges ( Figure 3). Most gauges only record one of the modeled tsunami events, although Fort Denison in Sydney Harbor records all five ( Figure 3). By using multiple sites with good quality data for each historic tsunami, we reduce the risk that site-specific factors limiting the model performance are mistakenly attributed to source/ propagation model biases (e.g., undiagnosed errors in the bathymetry or tide-gauge records). No gauge observations were rejected on the basis of disagreement with our tsunami model, to avoid biasing the results.
Tide-gauge data was obtained from a range of sources (see Acknowledgments). To detide the observations we first subtracted tidal predictions, which were either provided with the tide-gauge data, or obtained from TPXO7.2 (Egbert and Erofeeva, 2002). Next a high pass filter was used to remove the residual longperiod sea-level variations by applying a discrete Fourier transform and zeroing the amplitude of waves with period exceeding a threshold. At most sites the tides are semi-diurnal and the high pass filter threshold was 3 h. At our Western Australia sites the tides are diurnal and a longer threshold was used (6 h); the 3 h threshold led to non-stationary oscillations in the non-tsunami sea-level component at the time of the Sumatra 2004 tsunami, suggesting a longer threshold is appropriate. Irrespective this has little impact on the tsunami signal. In all cases the detided tsunami record might still be affected by other physical processes (e.g., seiching due to transient atmospheric forcings) or measurement errors (e.g., excess mechanical smoothing in the gauge, Satake et al., 1988); however it represents our best estimate of the tsunami signal.
Most gauges used herein have a sampling frequency of 1-min. For the Sumatra 2004 event, three of the four gauges used in Western Australia have 5-min sampling frequency (excepting Hillarys which has a 1-min frequency). For the 1960 Chile event only two gauges are available (Fort Denison and Cronulla), both of which were digitized at approximately 2-min sampling frequency from scans of the analogue tide-gauge records (the former by . Although in recent decades additional 15-min tide-gauge data is available at several sites, this was not used because it under-samples the tsunami making interpretation more difficult (Rabinovich et al., 2011). At Port Kembla the outer-harbour gauge is used; the inner-harbour gauge was not used based on advice from the data custodians (NSW Port Authority, personal communication 2020) and the issues noted in Allen and Greenslade (2016).

Tsunami Model
The tsunami is simulated from source to the tide-gauges using the open-source hydrodynamic model SWALS, which solves several variants of the shallow water equations in spherical or Cartesian coordinates and was used for the 2018 Australian Probabilistic Tsunami Hazard Assessment (Davies and Griffin, 2018). The source code includes a validation test-suite of more than 20 Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 598235 analytical, laboratory and field problems, including well-known tests such as those in NTHMP (2012) (excluding landslides) and two recent field-scale NTHMP problems (Macías et al., 2020). It uses structured grids with two-way nesting, and flux-correction to enforce conservation at nested grid boundaries. Different grids can employ different solvers simultaneously; this is used in the current study to test a range of reduced-physics solvers on the global grid (only) while solving the NSWE on refined grids. The global scale tsunami propagation is simulated on a 1 arcminute grid ( Figure 3A) using bathymetry derived from GEBCO 2014 (Weatherall et al., 2015) and GA250 (Whiteway, 2009) (details in Davies and Griffin, 2018). East-West periodic boundary conditions are used. The southern boundary (79°S) is covered by land. At the northern boundary (68°N) a reflective wall is imposed; while artificial this closes the model and facilitates energy conservation calculations (as SWALS does not track energy fluxes through boundaries). Physically this is reasonable given little tsunami energy will radiate through either Bering strait (which is small) or north Atlantic (which is very far from our tsunami sources). The reduced-physics solvers used on this grid are described in Sections 2.4.2-2.4.4.
Nearshore regions containing good-quality tide-gauge observations are simulated on refined grids (1/7 and 1/49 arcminutes, ( Figures 3C-E), which are linked with the global grid and each other via two-way nesting. Good quality nearshore elevation data is used on the highest-resolution grids, mostly derived from LIDAR and gap-filled using available single-beam bathymetric surveys and gridded elevation from prior studies (see Acknowledgments; Allen and Greenslade, 2016; . Breakwalls close to our tide-gauges were burned into the model's elevation to ensure they are represented irrespective of the grid resolution. On these nearshore grids the tsunami is simulated using the NSWE with Manningfriction (Section 2.4.1).
Results were compared with the default model at twelve sites with corresponding good-quality observations, as this is where the model results will be used. Waveforms were similar although not completely convergent, with tsunami maxima showing changes ranging from −12% to +21% with a median of 2% (typical examples in Figure 4). In general greater differences are anticipated between the models and data, so the default resolution was considered adequate for the current study.

Hydrodynamic Simplifications for Global Scale Tsunami Propagation
Several reduced-physics solvers are tested on the global grid, all based on the linear shallow water equations (LSWE) with alternative dissipation models: (1) Frictionless (Section 2.4.2) (2) Nonlinear Manning-friction (Section 2.4.3) (3) Constant linear-friction (Section 2.4.4) These represent computationally efficient alternatives to solving the full NSWE for global tsunami propagation. The global grid has a spatial resolution of 1-arc-minute and is modeled with reduced-physics approaches; (B) zoom near Australia which shows the three regions where refined grids are used and the NSWE are solved; (C-E) the refined grids around western and south-eastern Australia. Dashed boxes denote regions with 1/7-arc-minute spatial resolution (x260 m), and solid boxes denote regions with 1/49-arc-minute resolution (x 37 m) which contain the tide-gauges.
October 2020 | Volume 8 | Article 598235 The frictionless LSWE were tested because they are often used to simulate large-scale tsunami propagation (e.g., Choi et al., 2003;Burbidge et al., 2008;Thio, 2015;Li et al., 2016;Davies and Griffin, 2020) and arguably represent the simplest model of this kind. The Manning-friction approach was tested because it is closer to the full NSWE; however over timescales where dissipation is important, the nonlinearity will prevent a unitsource implementation. The constant linear-friction approach (Fine et al., 2013;Kulikov et al., 2014) was tested because it enables dissipation to be included in a linear framework, with all the associated efficiency benefits. Furthermore, although unitsources are not employed in this study, below it is shown that solutions of the constant linear-friction model can be well approximated with a simple transformation of solutions of the frictionless LSWE, which is convenient because the latter are available in several existing tsunami propagation databases (e.g., Thio, 2015;Li et al., 2016;Davies and Griffin, 2018).

Nonlinear Shallow Water Equations
The 2D NSWE with friction are (e.g., Baba et al., 2015;Behrens and Dias, 2015): is the 2D velocity vector (m/s), g is gravity (m/s 2 ), S f n 2 u u h −4/3 is the friction slope vector with Manning-friction coefficient n, and Ω ω(−vh, uh) is the Coriolis force with latitude dependent Coriolis parameter ω. If friction is neglected (S f 0) then for smooth flows Eq. 2 is energy conservative (Arakawa and Hsu, 1990): Here the depth-integrated energy density (e k + e p ) is integrated over a two-dimensional domain D with no inflow or outflow through boundaries, dA is an area element, E 0 is the constant total energy in D, e k is the depth-integrated kinetic energy density: with water density ρ (kg/m 3 ), and e p is the "depth-integrated available potential energy density": By definition the integral of e p in D is zero if the fluid mass is redistributed in D to make the free-surface constant (Lorenz, 1955); herein C 2 is a constant that satisfies this definition. If the vertical datum is the model's mean sea level, and there is no wetting and drying, then Eq. 5 simplifies to e p ρg 2 η 2 in wet regions and zero in dry regions. The latter form is common in the tsunami literature (Tang et al., 2012;Saito et al., 2014;Tolkova, 2014) but does not generalize to flows with wetting and drying or other vertical datums.
On refined grids the model herein solves Eq. 2 in fluxconservative form with a constant Manning-friction n 0.03. The refined-grid n value might be improved with site-specific tuning, but that was not attempted herein. A second-order accurate finite-volume scheme is used based on the hydrostatic-reconstruction approach of Audusse et al. (2004). The details are similar to the ANUGA software's DE1 flow algorithm (Davies and Roberts, 2015), although the latter uses a triangular mesh while the solver herein uses structured grids with two-way nesting.

Linear Shallow Water Equations Without Friction
In the deep ocean even large earthquake tsunamis are of small amplitude relative to ocean depths and velocities are slight compared to the gravity wave speed. Under these conditions Eq. 2 is well approximated with the frictionless LSWE: Here h 0 is the time-invariant depth below mean-sea-level, and the friction slope is zero but is included in Eq. 6 to facilitate extensions below. These equations are also energy conservative.
Herein the LSWE are solved numerically with the classical leap-frog scheme (Goto and Ogawa, 1997). Importantly this scheme has good numerical energy conservation; in our twelve simulations that used the frictionless LSWE on the global grid the numerically integrated energy (Eq. 3) varied within [−2.7%, 0.25%] over 60 h. Minor dissipation is expected in the full model because refined grids solve the NSWE with friction and employ a more dissipative finite-volume numerical scheme; minor energy increase can occur because the leap-frog scheme is not perfectly energy conservative. In comparison the literature suggests global tsunami propagation simulations often have much greater numerical dissipation; Tang et al. (2012) and Tolkova (2014) used several variants of the MOST and CLIFFS solvers to simulate the Tohoku tsunami globally without friction and found x 80-90% of the total energy was numerically dissipated in 24 h. Popinet (2011) reported substantial numerical dissipation when modeling the Sumatra 2004 tsunami from source with the Gerris finitevolume solver. Although the frictionless LSWE are energy conservative, physically some dissipation is expected. Considering the observed energy e-folding timescale of global tsunamis at late-times (x1 day) these frictionless LSWE simulations should under-estimate the physical dissipation, but the numerical scheme offers a good basis for testing dissipation models because it adds negligible numerical dissipation.

Linear Shallow Water Equations With Manning-Friction
To add dissipation to the global model it is natural to consider a Manning-friction closure in Eq. 6. For computational efficiency our approach exploits the approximation hxh 0 inherited from the LSWE: where n is the Manning-friction coefficient, which was set to 0.035 on the global grid (slightly larger than our nearshore Manning coefficient of 0.03). Real tsunamis likely dissipate due to a mixture of nonlinear bottom friction in shallow regions, and other deep-ocean dissipation mechanisms which are important for simulating tides (Egbert and Ray, 2000;Buijsman et al., 2015;Kleermaeker et al., 2019). Our global model poorly resolves most shallow nearshore areas, neglects nonlinear tidal interactions, and ignores other dissipation mechanisms, so the Manning model is likely a strong simplification of reality. In preliminary simulations we also tried a smaller global Manning coefficient (n 0.02) which did not perform as well; however further optimization of the global n may be possible.
Equation 8 is nonlinear so solutions cannot be constructed from unit-sources (at least not over timescales long enough for nonlinear dissipation to be important). However the term in large parenthesis is time-invariant which reduces the computational expense of our implementation. A semiimplicit discretization is used to include Eq. 8 in the leapfrog scheme, with q evaluated explicitly and q evaluated implicitly. The total energy dissipation over 60 h varied between 63 and 96% in our twelve simulations using this model on the global grid, with the lowest percentage dissipation corresponding to the smallest tsunami (2015 Chile) as expected with quadratic friction.
For global tsunami propagation this model is relatively close to the NSWE with Manning-friction (Eq. 2). Figure 5 compares solutions at some nearshore tide-gauges when using the above reduced-physics model on the global grid, vs. the full NSWE (the latter were solved with a leap-frog numerical scheme combined with an upwind treatment of nonlinear advection, for similarity with our reduced-physics solver, Goto and Ogawa, 1997;Liu et al., 1995). At our nearshore gauges of interest, which are inside the refined grids and thus simulated with the NSWE finite-volume scheme (Section 2.4.1), the difference between the two models is smaller than differences caused by grid refinement (compare Figures 4 and 5). However the full simulation ran three times faster when using the simpler model on the global grid (about 2 vs. 6 h for the 60 h simulation with all nested grids, using 8 CPU nodes of the Gadi supercomputer NCI, 2020), due purely to a factor-of-6.6 speedup in the global grid calculation.

Linear Shallow Water Equations With Constant Linear-Friction, and Approximation via Frictionless Solutions
To represent global-scale tsunami dissipation with a linear model, Fine et al. (2013) and Kulikov et al. (2014) appended constant linear-friction to the LSWE (Eq. 6) by replacing z zt → z zt + f in the momentum equations only, implying: where f is a constant linear drag coefficient. Fine et al. (2013) and Kulikov et al. (2014) show Eqs 6 and 9 cause the globally integrated energy to decay exponentially in time with e-folding timescale of 1/f , thus mimicking classical observations of tsunami energy time-decay at late times (Miller et al., 1962;Munk, 1963;van Dorn, 1984;van Dorn, 1987;Rabinovich et al., 2011;Rabinovich et al., 2013). On this basis they used f 1 × 10 − 5 which corresponds to an energy e-folding time of 27.7 h. Herein an implicit discretization is used to include Eq. 9 in the leap-frog scheme for Eq. 6. The total energy dissipation over 60 h was close to 88.5% for all twelve simulations using this model on the global grid, as expected theoretically. Linear-friction is attractive, if it works well enough in practice, because the model retains all the practical benefits of linearity for simulating global propagation (e.g., unit-source solutions are exact). The use of a constant f is not essential to preserve linearity; often linear friction is implemented by linearizing a quadratic drag model about a reference depth and velocity (e.g., Tolkova, 2015). However the use of a . To illustrate this point, Figure 6 compares a numerically derived solution to the constant linear-friction model Eqs 6 and 9 with an "approximate linear-friction" solution defined as: Here the superscript a denotes the approximate linear-friction solution, and the superscript 0 denotes the frictionless LSWE solution which is prescribed the same initial conditions as the solution of interest. Although Eq. 10 does not produce an exact solution of the constant linear-friction model, the error is negligible in practice ( Figure 6) for global-scale tsunamis and other waves having a period much smaller than the decay timescale 1/f . This is justified mathematically below. Thus unit-source databases based on the frictionless LSWE can be used, without modification, as a convenient means to derive solutions with linear-friction using Eq. 10. This makes the model very easy to use in practice, via existing frictionless tsunami scenario databases, and motivates us to test it herein.

Justification of the Approximate Linear-Friction Solution (Eq. 10)
It is standard to convert the frictionless LSWE to a single wave equation for the free surface by applying z zt to the mass equation and ∇· to the momentum equation in Eq. 6, and eliminating mixed derivatives of q. Using the superscript 0 to denote solutions of the frictionless LSWE, this leads to: where Ω 0 ω[−vh 0 , uh 0 ]. For solutions of the constant linearfriction model (denoted with superscript l) the same calculation yields an additional source-term on the RHS: For the approximate linear-friction solution (Eq. 10), the analogous equation follows by substituting [η 0 , uh 0 , vh 0 ] exp f 2 t [η a , uh a , vh a ] into Eq. 11: Equation 13 is the same as Eq. 12 except for the second RHS source term, which has magnitude xαf 2 /4 where α is the wave amplitude. If the wave has period T then the average absolute value of zη a zt is 4α/T. Thus the first RHS source-term has magnitude xf (4α/T), so is x1000 times larger than the second source term for a typical earthquake-generated tsunami in the deep ocean (period 30 min, f x10 − 5 ). In this circumstance Eq. 13 is a minor perturbation of Eq. 12 and their solutions will be arbitrarily close if fT/16 is sufficiently small, which explains the excellent agreement in practice ( Figure 6).

RESULTS
Section 3.1 provides examples to illustrate how the global dissipation models affect nearshore tsunami simulations in comparison with data. This gives context for a subsequent statistical analysis but necessarily represents a small sample of the results; figures comparing models and data at all sites are available via the repository (see Data Availability Statement). Section 3.2 uses statistical techniques to evaluate each model by comparing all simulations and observations simultaneously, focusing on biases in the modeled tsunami maxima and their relation to the simulation duration. These results motivate tests of several modified linear-friction models on the global grid (Section 3.3). Section 3.4 considers how the model performance varies among the source-inversions.

Fort Denison, Chile 1960, H19 Source Model
The 1960 Chile tsunami was observed widely on Australia's eastcoast where it induced widespread marine hazards and minor inundation (Beccari, 2009). Based on our simulations, the tsunami reached the eastern Australian mainland x 14 h after the earthquake via waves which propagated south of New Zealand. Further waves arrived via the ocean north of New Zealand, including due to prominent scattering of the leading wave around bathymetry near French Polynesia (x 11 h postearthquake, reaching eastern Australia x22 h post-earthquake). This led to a steady growth in the nearshore tsunami energy on Australia's east coast (Figure 7). At Fort Denison the largest waves occurred around 25-29 h post-earthquake, while they were slightly earlier (x22 − 25 h) at the Cronulla gauge which is 25 km south and closer to the open coast.
At Fort Denison the models with frictionless and Manningfriction global grid solvers give similar results prior to the observed maxima, with predicted waves slightly larger than observed (Figure 7). Linear-friction on the global grid produces slightly smaller initial waves. Constant linear-friction induces greater dissipation in the deep-ocean compared to Manning-friction which dissipates most energy in nearshore areas; at later times the cumulative effect of nearshore propagation leads to significant global energy dissipation with Manning-friction, but this has little influence on the early arrivals. All models well approximate the timing of the maximum wave but differ in its size (Figure 7). In this case linear-friction best represents the observed tsunami maxima, but at later times consistently under-predicts the tsunami size. Manning-friction better simulates the size of late arriving waves, while the frictionless global model leads to nearshore waves that are too large at late times.
The data has a phase-lag relative to all models (Figure 7). This is expected given our use of the LSWE on the global grid. The phase-lag should be better simulated by including additional terms in the global model related to loading, seawater density stratification, and self-gravitation (Watada et al., 2014;An and Liu, 2016), although this would be computationally expensive (Allgeyer and Cummins, 2014;Baba et al., 2017).

Hillarys, Sumatra 2004, F07 Source Model
The Sumatra 2004 tsunami induced significant marine hazards in Western Australia, including 35 ocean rescues, damages to boats and marinas, and minor inundation at a number of coastal towns (Anderson, 2015). It took about 6 h to reach Hillarys in Western Australia (Figure 8), which is the shortest travel time among the tsunamis considered herein. The fourth wave was the largest at Hillarys and occurred a few hours after the tsunami arrival, although waves around 21 h post-earthquake were only slightly smaller (Figure 8). Pattiaratchi and Wijeratne (2009) noted these later waves were the largest observed at several other sites in Western Australian, and attributed them to tsunami energy reaching Australia via reflections off distant Indian Ocean topography. At Hillarys the simulated leading waves are not strongly affected by the choice of global dissipation model (Figure 8). About 30 h post-earthquake the frictionless global model results in larger nearshore waves than simulations with friction, although initially it is not obvious that any model better agrees with data. However from about 36 h post-earthquake the frictionless global model consistently predicts larger nearshore waves than were observed, as seen in the previous Chile 1960 example, indicating a lack of global dissipation. At late times the models with global friction predict smaller waves and agree much better with data in the nearshore (Figure 8).
The data evidences some phase-lag relative to the model, as noted in the previous Chile 1960 example. For instance the modeled long-period wave around 20-23 h arrives slightly earlier than the observed wave (Figure 8), and a similar result is obtained with the P07 and L10 source models (see online repository). These phase-lags may be better modeled by accounting for loading, seawater density stratification, and self-gravitation (Baba et al., 2017). The neglect of wave dispersion is also likely to be significant for the Sumatra 2004 example, which shows much more short-period wave energy than the Chile 1960 example (e.g., compare Figures 7 and 8). Nondispersive shallow-water wave theory will over-estimate the celerity of shorter period waves, as previously noted in the context of the 2004 Sumatra tsunami (Kulikov, 2006).

Twofold Bay, Tohoku 2011, S13 Source Model
The 2011 Tohoku tsunami was widely observed in eastern Australia (Hinwood and Mclean, 2013). The initial waves reached Australia's east-coast via straits around eastern Papua New Guinea, the Solomon Islands and Vanuatu, while at later times much of the wave energy arrived via the ocean north of New Zealand (Hinwood and Mclean, 2013). Our models suggest the leading wave reached South America around 20-22 h postearthquake where it was prominently reflected into the Pacific Ocean, adding significantly to the late-time wave energy in the southwest Pacific Ocean. In eastern Australia the tide-gauges at Fort Denison and Twofold Bay exhibited late-time maxima (46-50 h post-earthquake, Figure 9) while gauges at Botany Bay and Port Kembla exhibited earlier maxima (18-20 h postearthquake) but nonetheless showed relatively large late-time waves that we attribute to reflections from the eastern Pacific. This is similar to reports from New Zealand where Borrero et al. (2015a) found some sites (but not all) experienced late-time tsunami maxima up to 2 days post-earthquake.
At Twofold Bay the first few waves (up to 20 h postearthquake) are similar in the models with frictionless and Manning-friction global grid solvers, and agree quite well with data ( Figure 9). Linear-friction produces smaller initial waves than the Manning model, as noted for above for the Chile 1960 tsunami, which is attributed to its greater dissipation in the deep ocean. At late times the frictionless global model produces overly large waves in the nearshore compared to data (as in the previous examples). In this case the Manning and linear-friction models predict late-time waves smaller than were observed. These results again highlight the sensitivity of later waves to the global dissipation model. Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 598235

Effect of Global Dissipation at Nearshore Gauges: Tsunami Maxima Statistics
To compare the models and observations at all gauges simultaneously the tsunami maxima (i.e. detided water-level maxima) is used following Allen and Greenslade (2016) and Adams and LeVeque (2017). The results above suggest that model biases may be related to the simulation duration, and to assess this both simulations and data are truncated to various time-intervals ( Figure 10): (1) 0-8 h after the tsunami arrival. The arrival time is defined, for each model and gauge separately, as the time that the modeled stage absolute value exceeds 5 × 10 − 4 of its maximum, and is typically 12 h post-earthquake herein (range of 5.8-16.0 h). (2) 0-36 h post-earthquake.
All but one observed time-series extends for the full 60 h simulation; the 1960 Cronulla tide-gauge record is truncated at 29 h post-earthquake and dropped from comparison with longer simulations. Statistics describing the model bias (G m ) and accuracy (|G| m ) are reported for each variation of the global model m ∈ (Frictionless, Manning, Linear) ( Figure 10). These emphasize the error as a fraction of the observation and deliberately weight each source inversion equally, even though they have different numbers of corresponding gauge observations, because biases in any one source-inversion will lead to correlated errors among its gauges (Section 3.4).
The bias statistic G m summarizes the relative model bias (e.g., G m −0.1 suggests 10% under-estimation is typical). It is calculated as: Here G m i gives the median relative error for model m and sourceinversion i (i 1 . . . 12) over all available tide-gauges j. The tidegauges have observed tsunami maxima d j and modeled tsunami maxima p m i,j . The superscript"**" is appended to G m values that show strong evidence of being significantly different from zero. It is applied when 10 or more of the 12 G m i values have the same sign. This criterion is heuristic but is motivated by the following argument: if the model m and source-inversions i have little bias on average then any G m i has an equal chance of being positive or negative. If all G m i are independent then the signs of the G m i behave like a binomial random variable, and with 12 source models the resulting binomial distribution (parameters p 0.5, size 12) implies a 96% chance that fewer than 10 have the same sign.
The accuracy statistic |G| m is similar to G m , except the absolute value of the relative model error is used: Thus |G| m estimates the typical magnitude of the model error as a fraction of the data, irrespective of its sign. Figure 10 highlights the interaction of the simulation duration and the model bias, most prominently for the frictionless global model. If simulations are restricted to 0-8 h following tsunami Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 598235 arrival, this model exhibits small bias and a typical accuracy around 24% (top-left panel of Figure 10). For a 36 h simulation there is increasing positive bias (x24%), which becomes pronounced in the 60 h simulation (x43%). The frictionless model bias on the last day of 60 h simulation is very strong (x82%), emphasizing that it consistently overestimates the size of late-arriving waves (Section 3.1).
Manning-friction in the global model leads to much better agreement with observed tsunami-maxima for long simulations (Figure 10). For short simulations (0-8 h post-arrival) it FIGURE 10 | Modelled-vs-observed tsunami maxima for different global model types (frictionless, Manning-friction, linear-friction), and different temporal subsets of the simulation (0-8 h after tsunami arrival; 0-36 h after earthquake; 0-60 h after earthquake, 36-60 h after earthquake). Diagonal solid line is y x; diagonal dashed lines are y 1.5x and y x/1.5. Where the G m values are followed by **, there were 10 or more source-inversions having G m i values of the same sign.
Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 598235 13 performs similarly to the frictionless model, suggesting friction has limited influence on early arriving waves. But over long simulations the cumulative effect of nearshore dissipation at the global scale reduces the nearshore wave heights, in a manner broadly consistent with observations ( Figure 10).
Linear-friction in the global model induces a moderate negative bias in the early (0-8 h post-arrival) nearshore tsunami maxima (x − 28%**, Figure 10). The constant linearfriction formulation induces greater deep-ocean dissipation than Manning-friction, and the results here suggest this dissipation is too large. However the bias weakens for longer simulations, and for the final 36-60 h post-earthquake the model performs as well as Manning-friction ( Figure 10). Compared with the frictionless model, which shares the advantages of linearity, the constant linear-friction model is clearly superior at late times but inferior at short times ( Figure 10).
To gauge the model accuracy it is useful to cross-reference the G m and |G| m statistics in Figure 10 with previous studies which compared modeled and observed tsunami maxima for multiple source-inversions. Allen and Greenslade (2016) modeled 9 historic tsunamis at Port Kembla outer-harbour with the MOST model, using source inversions based on the T2 database (Greenslade et al., 2011). From data in table 3 of Allen and Greenslade (2016) we computed G m −0.1 and |G| m 0.36. Their results at another inner-harbour site were less accurate but are not reported herein because the tidegauge accuracy is doubtful (NSW Port Authority, personal communication 2020). Adams and LeVeque (2017) studied five source-inversions with 4-6 gauges each using two models (GEOCLAW and MOST). Using results tabulated in Figure 3 of Adams and LeVeque (2017) we compute G m −0.11, − 0.01 and |G| m 0.375, 0.213, respectively. Given the small number of source-inversions and gauges used in all studies, and their different methodologies, differences between these statistics are probably not meaningful. But they suggest model errors comparable to those obtained with our methodology using global-scale Manning-friction (at all times); using the globalscale frictionless model (at early times); and using global-scale linear-friction (at late times).

Alternative Linear-Friction Models With Reduced Bias for Earlier Waves
Because linear models have significant practical benefits (e.g., unit sources), herein two variations of the linear-friction model are tested. Both aim to reduce downward bias in nearshore waves 0-8 h after arrival, as compared with the constant linear-friction model, while retaining the benefits of friction for longer simulations. In addition, both models retain the convenient property that their solutions can be approximated from solutions of the frictionless model: (1) Reduced-linear-friction.
This model uses f 1/(36 × 3600)x7.71 × 10 − 6 , corresponding to an artificially long tsunami energy e-folding timescale of 36 h (vs. 27.7 h in the original model). The 36 h decay-timescale will reduce energy loss but is a-priori expected to be too long, noting Rabinovich et al. (2013)  (2) Delayed-linear-friction. This model is frictionless for the first 12 h, and subsequently uses f 1 × 10 − 5 as for the original linear-friction model. The heuristic motivation for an initial frictionless period is that early waves reach Australia predominantly via the deep-ocean, with few nearshore interactions, so the cumulative effect of bottom friction should initially be small. The 12 h time-period matches the median tsunami arrival time for the events studied herein (Section 4 will consider alternative estimates of this time-period). In contrast to reduced-linear friction, this model's late-time energy-decay timescale remains broadly consistent with observations (Rabinovich et al., 2013). Solutions of this model can still be well approximated by frictionless solutions via straightforward modification of Eq. 10: where t 12 gives the number of seconds in 12 h. Figure 11 compares the performance of these models and the original linear-friction model with data. The reduced-linear-friction model continues to show downward bias for waves 0-8 h after arrival, even though its friction coefficient is a-priori too small ( Figure 11). If constant linear-friction were a good approximation of dissipation for early-arriving waves, then the a-priori small friction coefficient should lead to overestimation of tsunami maxima, but this does not occur. That indicates weaknesses in the constant linearfriction parameterization of tsunami dissipation; friction is still too large in the deep ocean. However the model performs reasonably well for longer simulations (Figure 11).
In contrast, the delayed-linear-friction model performs better in the 0-8 h post-arrival period than any other linear model with friction ( Figure 11). Furthermore, over longer simulations it continues to exhibit relatively low-bias in the tsunami maxima, with results comparable to the Manning-friction model ( Figure 10). Among the linear models tested herein, this appears the most promising overall for simulating globalscale tsunami propagation.

Source-Inversion Effects on Modeled Tsunami Maxima
The previous section highlights that error in the modeled tsunami maxima varies from gauge-to-gauge and also depends on the simulation duration. However there is also a significant component related to the source-inversion itself ( Figure 12). For simplicity Figure 12 only depicts results using the Manningfriction and delayed-linear-friction global models, which exhibited the least overall bias above, although comparable variations between source-inversions exist for the other models.
For a given historical tsunami, inversions with a larger model/ observed ratio in Figure 12 tend to have greater available potential energy in their ocean-displacement (reported in Figure 2). For example consider the Chile 2010 tsunami; the L11 inversion produces higher model/obs ratios than the F13 inversion, while the former also exhibits greater available potential energy (Figures 2, 12). For the Chile 2015 tsunami the W17 inversion exhibits greater model/obs than the R16 inversion, and has greater initial available potential energy (Figures 2, 12). Other cases behave similarly with one FIGURE 11 | Modelled-vs-observed tsunami maxima for different global model types (linear-friction, linear-friction with a lower friction factor, linear-friction with delayed onset of 12 h), and different temporal subsets of the simulation (0-8 h after tsunami arrival; 0-36 h after earthquake; 0-60 h after earthquake, 36-60 h after earthquake). Diagonal solid line is y x; diagonal dashed lines are y 1.5x and y x/1.5. The LHS column is the same as in Figure 10.
Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 598235 15 exception; for Sumatra 2004 the P07 inversion has slightly smaller model/observed than the F07 inversion ( Figure 12), despite having a slightly larger available potential energy ( Figure 2). This may reflect that compared to F07, the P07 inversion has greater co-seismic displacement toward the north of the rupture area (Figure 2), which may have less effect on Australia than deformation in the southeast due to tsunami dynamics in the Bay of Bengal. In any case this result is consistent with the suggestion of Titov et al. (2016) that farfield tsunamis are sensitive to the location and energy of the initial ocean-surface displacement, with other source details having a secondary role.

DISCUSSION
Linear models provide a very convenient basis for global-scale tsunami scenario databases, and these are often implemented without friction (Burbidge et al., 2008;Li et al., 2016;Davies and Griffin, 2020). Our results suggest this approach is adequate for simulating tsunami maxima up to 8 h following tsunami arrival, when combined with nearshore models based on the NSWE with Manning-friction. Herein 8 h post-arrival corresponds to a median of 20 h post-earthquake; for hazard applications this duration would typically include the most significant waves for the most hazardous scenarios, where sites of interest are well exposed to the tsunami's leading wave. However, tsunamis may remain dangerous for a considerably longer duration. Over longer simulations the lack of global model dissipation leads to slow divergence of the frictionless model wave heights, as compared to both data and models with friction. It is difficult to estimate a "threshold" duration after which the frictionless model biases become important. In comparison to observed tsunami maxima, our results for 36 h simulations suggest bias x24% in the frictionless model, but considering variations among the sourceinversions this is not strong enough to be conclusive. However the biases are very clear in frictionless simulations of 60 h duration, especially if the model-data comparison is restricted to final 36-60 h post-earthquake.
The late-time model biases are largely removed using nonlinear Manning-friction in the global propagation model with n 0.035. For simulations of 24 h duration the resulting tsunami maxima are often x 10% smaller than with the frictionless model, but the difference grows for longer simulations (top row of Figure 13). Because Manningfriction is nonlinear this approach cannot be applied using unit-sources (at least not over timescales long enough for nonlinear-dissipation to be important), but it is a good option if one can simulate the tsunami from source. From a physical perspective this model remains over-simplified; it neglects tides which will interact nonlinearly with the tsunami, and furthermore, the tidal literature suggests about one-third of tidal-dissipation occurs in the deep ocean due to mechanisms that are not represented with shallow-water bottom friction (Egbert and Ray, 2000;Buijsman et al., 2015). However we focus on pragmatic methodologies for applications, and by this standard Manning-friction is a good choice for late-time global-propagation simulation, assuming it is viable to model tsunami propagation from source.
In some instances unit-source based treatments of global tsunami propagation are essential; for example if many tsunami scenarios must be simulated for a small nearshore area, global propagation modeling may be computationally prohibitive. In such circumstances the delayed-linear friction model is worth considering, particularly for long simulation durations where biases in the frictionless model are anticipated. While not motivated from hydrodynamic theory, this model addresses the early-time bias of constant linear-FIGURE 12 | Per-inversions boxplots of ratio between modeled and observed tsunami maxima at all gauges, for simulations of 60 h, and two model types. The integers above the boxes count the tide-gauge observations. Dotted horizontal lines are y 1.5 and y 1/1.5. Definition of the boxplots: The box defines the interquartile range with the central line defining the median. The "whiskers" extend from the box for at most 1.5 inter-quartile ranges, or up to the range of the data. Points not covered by whiskers are plotted directly. The Chile 1960 sources only have one gauge record that covers the 0-60 h simulation (Figure 7), so the boxplot reduces to a horizontal line.
Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 598235 friction, which was found to underestimate observed tsunami maxima 0-8 h post-arrival. A key result of this study is that the model can be implemented, to high-accuracy, via a trivial transformation of frictionless solutions (Eq. 16). Thus the model is very easy to implement with existing unit-source databases based on the frictionless LSWE (e.g., Thio, 2015;Davies and Griffin, 2020).
Herein delayed-linear-friction was applied with a 12 h frictionless period (t 12 ), and it is natural to ask if this duration is optimal. Figure 13 (bottom row) shows that while delayedlinear-friction generally predicts tsunami maxima similar to the Manning-friction model, there is a negative correlation between the modeled tsunami arrival time and the ratio of the two models ( Figure 13). This result is expected if the Manning-friction model is better approximated using the gauge-specific arrival time t a to define the frictionless period in Eq. 16 rather than t 12 , i.e.: One may estimate the relative change r in the tsunami maxima that would result from this approach, vs. the use of t 12 herein, by assuming the modeled maxima responds linearly to the incoming wave. In this case r is equal to the ratio of Eqs 16 and 17: where t M indicates the timing of the offshore waves controlling the tsunami maxima. The inequality constraints in Eq. 18 should hold for almost all of our delayed-linear-friction modeled gaugerecords; only 2/68 modeled time-series in Figure 13 have tsunami maxima arriving before 14 h post-earthquake, or less than 3 h post-tsunami-arrival; only one observed gauge maxima occurs before 18 h post-earthquake (Hillarys for Sumatra 2004; Figure 8). The fact that 1/r reasonably approximates the tsunami maxima ratio in Figure 13 suggests that, for sitespecific studies, the Manning-friction results may be better approximated with delayed-linear-friction using a "local" initial frictionless period defined by t a . This "local" approach was not implemented herein because a wide range of sites were modeled simultaneously, each having their own t a (Figure 3). However it would be straightforward to apply in site-specific studies, and merits further testing.
FIGURE 13 | Ratios of modeled tsunami maxima, at nearshore sites with tide-gauge data, using different offshore friction treatments: Top row frictionless/Manning; Bottom row Delayed-linear-friction/Manning-friction. In the bottom row the dashed line is y 1/exp(−f (t 12 − t a )/2) where t a is the modeled gauge-specific tsunami arrival time and t 12 is the number of seconds in 12 h; this is a linear approximation of the proportionate reduction in the delayed-linear-friction tsunami maxima if the initial frictionless period was t a instead of t 12 (under various assumptions appropriate for our sites, as discussed in the text).
Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 598235 Finally we emphasize limitations of this study and directions for further work. All results are based on globalscale propagation models with little numerical dissipation. Clearly these techniques should not be applied to models where numerical dissipation is already comparable to or greater than the required physical dissipation, as this could exaggerate any under-estimation of late arriving waves. Our models neglect various physical processes which can affect tsunami propagation (dispersion, loading, self-gravitation, density stratification), and we make no attempt to resolve coastal inundation outside our areas of interest in Australia although this may affect remote wave reflections; future work should consider the influence of these approximations on modeled waves at the coast. Following previous work (Allen and Greenslade, 2016;Adams and LeVeque, 2017) we used the tsunami-maxima statistic to quantify model biases, and while this is useful, there is much more information available in the full time-series that may be extracted using other statistics and provide new insights into the model performance. The tests herein consider only five historic tsunamis using sources derived from 12 finite-fault inversions, with a total of 28 tide-gauge records. The testing should be extended to consider more historic events and sites, and other types of data (e.g., inundation footprints and depths; tsunami currents). Of particular interest is the global-record of historical tsunami energy-decay timescales (e.g., van Dorn, 1987;Rabinovich et al., 2011Rabinovich et al., , 2013; further insights into dissipation models would likely be obtained by simulating these observations with alternative models. The tests should also be extended to consider tsunami scenarios derived for hazard assessment (Davies, 2019), including characterization of the nearshore performance of random hazard scenarios.

CONCLUSION
When modeling far-field tsunamis in the nearshore it is often convenient to nest high-resolution nonlinear shallow water models (which cover coastal regions of interest) within reduced-physics global-scale models (which efficiently simulate tsunami propagation). In this context we evaluated the performance of several reduced-physics global propagation models which combine the LSWE with alternative treatments of friction. Nearshore tsunamis were simulated with a nested nonlinear model and compared with coastal tide-gauge observations in Australia. Tsunami initial conditions were derived from published source-models. Our results suggest the commonly used frictionless global-scale model is adequate for simulating far-field coastal tsunamis for 0-8 h following arrival. However for sufficiently long simulations the frictionless model overestimates coastal wave heights because it is mathematically energy conservative (which is well approximated with our numerical methods) and thus does not represent global-scale tsunami dissipation. In our simulations this bias was clear from comparison with data 36-60 h post-earthquake; it is difficult to define a precise time at which the bias becomes significant. Better estimates of latetime wave heights can be obtained by appending Manningfriction to the global model, although this renders the model nonlinear. For unit-source based implementations it is preferable work with linear models, and so variants of linear-friction were tested; these are not derived from hydrodynamic theory yet can mimic observed tsunami energy decay rates. Among these, delayed-linear-friction most accurately simulated tsunami maxima at nearshore gauges. Solutions of this model can be conveniently derived by transforming frictionless LSWE solutions with Eq. 16, making it easy to implement via existing unit-source databases that are derived with the frictionless LSWE (e.g., Davies and Griffin, 2020). These results may facilitate improved simulation of late-time tsunami wave heights for hazard and early warning applications in the far-field.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/ GeoscienceAustralia/ptha/tree/master/misc/nearshore_testing_ 2020.

AUTHOR CONTRIBUTIONS
GD led the study design and implementation and wrote the first draft manuscript, with guidance from FR and SL. FR provided code and data to construct the sources R14 and R16. All authors contributed to the final manuscript.