Omori’s law states that the rate of aftershocks decays as a function of inverse time. There are multiple physical explanations that we reduce into a nonlinear mixed effects relation of three terms: (1) a Rate/State expression that can account for static/dynamic and viscoelastic triggering caused directly by the mainshock, (2) a fluid diffusion triggering term, and (3) a randomized secondary triggering (cascade) term. We fit free physical-model parameters to an observed aftershock sequence through two nonlinear regression methods to find the relative contributions of physics-based models in an observed aftershock sequence. Results from both methods show that Rate/State models overpredict aftershock rates by ∼0–30%. Secondary aftershocks cause a net negative contribution (seismicity rate reduction that corrects overprediction by other terms) ranging between ∼0 and 30%. All regression solutions yield negative secondary triggering contributions without being guided to do so. A physical explanation for this is that aftershock occurrence relieves stress from the crust, ultimately causing the sequence to extinguish itself. Fluid diffusion triggering contributions range from ∼0 to 20%. Diffusion processes are observed to be shorter in time than the full duration of an aftershock sequence and they are also spatially limited, diminishing their influence. Our results apply to an aftershock decay curve from the 2016 Central Apennines earthquake sequence, meaning that our specific results may not be general. Our primary conclusion is that any one physical model cannot alone fit the observed sequence as well as the combination of three we investigated.

We use combinatorial optimization to find the optimal spatial distribution of random samples of earthquakes (≥6.5) that minimize the misfit in target slip rates for all faults in the northeast Caribbean, and we derive magnitude–frequency relationships with uncertainties for these faults. Slip rates for many faults are derived from geodetic block models, not direct measurements, because of their underwater locations. Predicted recurrence rates for earthquakes on the East Hispaniola and Puerto Rico trench faults are 220–450 yr for moment magnitude (M) 7 and 3000–5000 yr for M 8, with the maximum feasible magnitude of M 8.2. The most frequent earthquakes with magnitudes ≥7.0 are predicted on the large upper plate strike-slip faults, Enriquillo (EF), and Septentrional faults, commensurate with the historical record. Calais et al. (2023) suggested that shortening in western Hispaniola is accommodated on the offshore Jérémie and onshore Malpasse faults north and south of EF, instead of on terrestrial faults in the western Hispaniola and EF. Because of our system-modeling approach, such a configuration predicts less frequent earthquakes on EF and western Hispaniola and Muertos convergent zones. Recurrence times of a few hundred years for M 6.7 earthquakes are predicted on the submerged faults in Mona Passage, and infrequent M > 7 earthquakes are predicted on the Virgin Islands faults.

As probabilistic tsunami hazard analysis (PTHA) focuses more on assessments for localized, populous regions, techniques are needed to identify a subsample of representative earthquake ruptures to make the computational requirements for producing high-resolution hazard maps tractable. Moreover, the greatest epistemic uncertainty in seismic PTHA is related to source characterization, which is often poorly defined and subjective. We address these two salient issues by applying streamlined earthquake rupture forecasts (ERFs), based on combinatorial optimization methods, to an unsupervised machine learning workflow for identifying representative ruptures. ERFs determine the optimal distribution of a millennia-scale sample of earthquakes by inverting the observed slip rate on major faults. We use two previously developed combinatorial optimization ERFs, integer programming and greedy sequential, to produce the optimal location of ruptures with seismic moments sampled from a regional Gutenberg-Richter magnitude-frequency distribution. These ruptures in turn are used to calculate peak nearshore tsunami amplitude, using computationally efficient tsunami Green's functions. An unsupervised machine learning workflow is then used to identify a small subsample of the earthquakes input to ERFs for onshore PTHA analysis. We eliminate epistemic uncertainty related to source distribution under traditional PTHA analysis; in its place, a quantifiable, less subjective and generally smaller uncertainty related to the input to ERFs is included. The Nankai subduction zone is used as a test case, where previous ERFs have been conducted. Results indicate that the locations of representative earthquakes are sensitive to choice of magnitude-area relation and to whether a minimum cumulative stress objective is imposed on the fault. In general, incorporating ERFs into PTHA provide a physically self-consistent method to incorporate fault slip information in determining representative earthquakes for onshore PTHA, eliminating a major source of epistemic uncertainty.

ABSTRACT We determine optimal on-fault earthquake spatial distributions using a combinatorial method that minimizes the long-term cumulative stress resolved on the fault. An integer-programming framework was previously developed to determine the optimal arrangement of a millennia-scale earthquake sample that minimizes the misfit to a target slip rate determined from geodetic data. The resulting cumulative stress from just slip-rate optimization, however, can greatly exceed fault strength estimates. Therefore, we add an objective function that minimizes cumulative stress and broad stress constraints to limit the solution space. We find that there is a trade-off in the two objectives: minimizing the cumulative stress on a fault within fault strength limits concentrates earthquakes in specific areas of the fault and results in excursions from the target slip rate. Both slip-rate and stress objectives can be combined in either a weighted or lexicographic (hierarchical) method. Using a combination of objectives, we demonstrate that a Gutenberg–Richter sample of earthquakes can be arranged on a constant slip-rate finite fault with minimal stress and slip-rate residuals. We apply this method to determine the optimal arrangement of earthquakes on the variable slip-rate Nankai megathrust over 5000 yr. The sharp decrease in slip rate at the Tokai section of the fault results in surplus cumulative stress under all scenarios. Using stress optimization alone restricts this stress surplus to the northeast end of the fault at the expense of decreasing the slip rate away from the target slip rate at the southwest end of the fault. A combination of both slip-rate and stress objectives provides an adequate fit to the data, although alternate model formulations for the fault are needed at the Tokai section to explain persistent excess cumulative stress. In general, incorporating stress objectives and constraints into the integer-programming framework adds an important aspect of fault physics to the resulting earthquake rupture forecasts.

Evaluation of tsunami disaster risk for a coastal region requires reliable estimation of tsunami hazard, for example, wave amplitude close to the shore. Observed tsunami data are scarce and have poor spatial coverage, and for this reason probabilistic tsunami hazard analysis (PTHA) traditionally relies on numerical simulation of "synthetic" tsunami generation and propagation toward the coast. Such an approach has been extensively studied in the past and it is widely recognized as an important disaster-risk mitigation tool. PTHA can not only provide less uncertain and spatially coherent hazard estimates in comparison with classical empirical data analysis which is restricted at the tide gauge stations, but also local inundation information. In this paper, we explore a purely statistical alternative to traditional PTHA for evaluation of tsunami amplitude hazard. Here, we use tide gauge measurements of tsunami amplitude along the western United States, specifically California and Oregon, and develop a spatial Bayesian hierarchical model (BHM) to assess tsunami hazard from far-field earthquake sources at various recurrence intervals. The configuration of our model incorporates latent Gaussian fields that utilize information on the distance between tide gauges as well as on the continental shelf width, that is, a covariate linked to potential dissipative effects on wave energy as the tsunami travels over shallow water. Through our BHM, we produce spatially continuous probabilistic maps of far-field tsunami hazard which can aid comprehensive tsunami disaster risk reduction and management. Devastating tsunamis are rare, but their consequences can be destructive for people who live close to the coast and their livelihoods. Assessing tsunami hazard, for example, nearshore tsunami height, typically requires us to resort to physics-based models since historical data from individual monitoring stations are scarce. In this work, we develop a spatial statistical model which can pool empirical amplitude data and capture dependence between tide gauges to allow for a more robust tsunami hazard analysis than analyzing station data in isolation. Our model can serve as a reliable benchmark and supplement to traditional probabilistic tsunami hazard analysis. A spatial Bayesian hierarchical model is proposed for tsunami amplitude data from far-field earthquake sources in the US West CoastLatent Gaussian processes are utilized to capture dependence between tide gauges through pair-wise distances and the continental shelf widthThe model allows for obtaining spatially continuous far-field tsunami hazard estimates at low probabilities of exceedance

We study stress-loading mechanisms for the California faults used in rupture forecasts. Stress accumulation drives earthquakes, and that accumulation mechanism governs recurrence. Most moment release in California occurs because of relative motion between the Pacific plate and the Sierra Nevada block; we calculate relative motion directions at fault centers and compare with fault displacement directions. Dot products between these vectors reveal that some displacement directions are poorly aligned with plate motions. We displace a 3D finite-element model according to relative motions and resolve stress tensors onto defined fault surfaces, which reveal that poorly aligned faults receive no tectonic loading. Because these faults are known to be active, we search for other loading mechanisms. We find that nearly all faults with no tectonic loading show increase in stress caused by slip on the San Andreas fault, according to an elastic dislocation model. Globally, faults that receive a sudden stress change respond with triggered earthquakes that obey an Omori law rate decay with time. We therefore term this class of faults as "aftershock faults." These faults release 4% of the moment release in California, have 0.1%-5% probability of M 6.7 earthquakes in 30 yr, and have a 0.001%-1% 30 yr M 7.7 probability range.

We use amplitude ratios from narrowband-filtered earthquake seismograms to measure variations of seismic attenuation over time, providing unique insights into the dynamic state of stress in the Earth’s crust at depth. Our dataset from earthquakes of the 2016–2017 Central Apennines sequence allows us to obtain high-resolution time histories of seismic attenuation (frequency band: 0.5–30 Hz) characterized by strong earthquake dilatation-induced fluctuations at seismogenic depths, caused by the cumulative elastic stress drop after the sequence, as well as damage-induced ones at shallow depths caused by energetic surface waves. Cumulative stress drop causes negative dilatation, reduced permeability, and seismic attenuation, whereas strong-motion surface waves produce an increase in crack density, and so in permeability and seismic attenuation. In the aftermath of the main shocks of the sequence, we show that the M ≥ 3.5 earthquake occurrence vs. time and distance is consistent with fluid diffusion: diffusion signatures are associated with changes in seismic attenuation during the first days of the Amatrice, Visso-Norcia, and Capitignano sub-sequences. We hypothesize that coseismic permeability changes create fluid diffusion pathways that are at least partly responsible for triggering multi-mainshock seismic sequences. Here we show that anelastic seismic attenuation fluctuates coherently with our hypothesis.

A critical component of seismic hazard analysis is understanding the frequency and spatial distribution of earthquakes with different magnitudes on nearby faults. A framework for determining the optimal spatial distribution of earthquakes on a complex fault system is developed using combinatorial optimization methods. Input to the framework is a millennia-scale sample of earthquakes taken from a regional Gutenberg-Richter (G-R) relation. We then determine the optimal spatial arrangement of each earthquake in the fault system according to an objective function and constraints. Our previously published results focus on minimizing the total misfit in slip rates as the objective function; constraints were maximum and minimum slip rate values that incorporate uncertainty in slip-rate values for each fault. Both global and local combinatorial optimization methods have been developed to solve these problems: integer programming and the greedy sequential algorithm, respectively. Resulting on-fault magnitude distributions cannot be simply classified as being either purely characteristic or G-R. For example, faults may exhibit multiple “characteristic” magnitudes or a power-law distribution of magnitudes over a restricted range. Current research involves adapting the general combinatorial framework to include other and multiple objective functions, including minimizing the variation in accumulated stress over millennia. The framework can also accommodate branching and step-over connections for the slip-rate objective, while current research is underway to include interaction stress loading among the different faults in the fault system for stress-based objectives. Results from these methods are valuable for verifying the assumed magnitude-frequency distributions for faults in probabilistic seismic and tsunami hazard analyses.

On‐fault earthquake magnitude distributions are calculated for northern Caribbean faults using estimates of fault slip and regional seismicity parameters. Integer programming, a combinatorial optimization method, is used to determine the optimal spatial arrangement of earthquakes sampled from a truncated Gutenberg‐Richter distribution that minimizes the global misfit in slip rates on a complex fault system. Slip rates and their uncertainty on major faults are derived from a previously published GPS block model for the region, with fault traces determined from offshore geophysical mapping and previously published onshore studies. The optimal spatial arrangement of the sampled earthquakes is compared with the 500‐year history of earthquake observations. Rupture segmentation of the subduction interface along the Hispaniola‐Puerto Rico Trench (PRT) fault and seismic coupling on the PRT fault appear to exert the primary control over this spatial arrangement. Introducing a rupture barrier for the Hispaniola‐PRT fault northwest of Mona Passage, based on geophysical and seismicity observations, and assigning a low slip rate of 2 mm/yr on the PRT fault are most consistent with historical earthquakes in the region. The addition of low slip‐rate secondary faults as well as segmentation of the Hispaniola and Septentrional strike‐slip fault improves the consistency with historical seismicity. An important observation from the modeling is that varying the slip rate on the PRT fault and different segmentation scenarios result in significant changes to the optimal magnitude distribution on faults farther away. In general, optimal on‐fault magnitude distributions are more complex and inter‐dependent than is typically assumed in probabilistic seismic hazard analysis and probabilistic tsunami hazard analysis.

The NEAM Tsunami Hazard Model 2018 (NEAMTHM18) is a probabilistic hazard model for tsunamis generated by earthquakes. It covers the coastlines of the North-eastern Atlantic, the Mediterranean, and connected seas (NEAM). NEAMTHM18 was designed as a three-phase project. The first two phases were dedicated to the model development and hazard calculations, following a formalized decision-making process based on a multiple-expert protocol. The third phase was dedicated to documentation and dissemination. The hazard assessment workflow was structured in Steps and Levels. There are four Steps: Step-1) probabilistic earthquake model; Step-2) tsunami generation and modeling in deep water; Step-3) shoaling and inundation; Step-4) hazard aggregation and uncertainty quantification. Each Step includes a different number of Levels. Level-0 always describes the input data; the other Levels describe the intermediate results needed to proceed from one Step to another. Alternative datasets and models were considered in the implementation. The epistemic hazard uncertainty was quantified through an ensemble modeling technique accounting for alternative models’ weights and yielding a distribution of hazard curves represented by the mean and various percentiles. Hazard curves were calculated at 2,343 Points of Interest (POI) distributed at an average spacing of ∼20 km. Precalculated probability maps for five maximum inundation heights (MIH) and hazard intensity maps for five average return periods (ARP) were produced from hazard curves. In the entire NEAM Region, MIHs of several meters are rare but not impossible. Considering a 2% probability of exceedance in 50 years (ARP≈2,475 years), the POIs with MIH >5 m are fewer than 1% and are all in the Mediterranean on Libya, Egypt, Cyprus, and Greece coasts. In the North-East Atlantic, POIs with MIH >3 m are on the coasts of Mauritania and Gulf of Cadiz. Overall, 30% of the POIs have MIH >1 m. NEAMTHM18 results and documentation are available through the TSUMAPS-NEAM project website (http://www.tsumaps-neam.eu/), featuring an interactive web mapper. Although the NEAMTHM18 cannot substitute in-depth analyses at local scales, it represents the first action to start local and more detailed hazard and risk assessments and contributes to designing evacuation maps for tsunami early warning.

Because of their inaccessibility, submarine landslides are typically studied individually and at great effort and expense to provide knowledge of the specific site conditions where these landslides occur. Statistical analysis of submarine landslide scars can offer generalized perspectives on the processes that initiate submarine landslides and can help toward hazard assessment in areas that have not been studied in detail. The following review discusses more than a decade of development of statistical approaches to studying submarine landslides. Landslides were previously viewed together with other natural hazards, such as earthquakes and fires, as a phenomenon whose size distribution obeys an inverse power law. Inverse power law distributions are the result of self-organized avalanche processes, in which the final hazard size cannot be predicted at the onset of the disturbance. We find that volume and area distributions of submarine landslides along the U.S. Atlantic continental slope and along nine other margins worldwide do not follow an inverse power law. Rigorous statistical tests of several different probability distribution models indicate that the lognormal model is most appropriate for these siliciclastic environments. Lognormal distributions can be simulated by assuming that the area of slope failure depends on earthquake magnitude, in other words, failure occurs simultaneously over the area affected by horizontal ground shaking and does not cascade from nucleating sources. Therefore, the maximum landslide size can be predicted from the earthquake magnitude and the distance from the rupturing fault. Moreover, earthquakes <~M4.5 cannot generate significant submarine landslides. We further demonstrate that empirical, offshore landslide hazard curves can be developed from these lognormal landslide size distributions, if the duration of mapped landslide activity is known. In addition to hazard estimation, scaling relationships can yield insights on the physical processes associated with landslide failure. For example, the log-log relationship between volume and area of landslide scars in siliciclastic margins is observed to be almost linear implying that most landslides are translational. Carbonate margins, in contrast, show a power-law distribution of scar volumes and their volume to area relationship is ~1.3. These results suggest that landslides in carbonate margins are governed by the random distributions of existing fissures, and they act like rock falls on land. Although earthquakes are the principal trigger of submarine landslides, the effects of earthquake frequency on slope stability can be counterintuitive. The average size of landslide scars decreases non-linearly with increasing frequency of earthquakes and increases with increasing sedimentation rate. The effect is interpreted as evidence for densification and shear strength increase of margin sediment, induced by repeated seismic shaking.

Abstract A new global optimization method is used to determine the distribution of earthquakes on a complex, connected fault system. The method, integer programming, has been advanced in the field of operations research but has not been widely applied to geophysical problems until recently. In this application, we determine the optimal distribution of earthquakes on mapped faults to minimize the global misfit in slip rates for multifault ruptures. Integer programming solves for a decision vector composed of every possible location that a sample of earthquakes can occur on every fault, subject to slip rate uncertainty constraints. Step over connections are straightforward to include, whereas branching fault connections are not. To include branching ruptures, we distinguish between individual multifault rupture paths, as opposed to formulating the integer programming problem based on individual faults as in previous studies. The new method is applied to the complex fault system in the San Francisco Bay Area as a case study. Results from the integer programming method are compared to those from a local optimization method, termed the greedy‐sequential method. Several experiments using these two methods indicate that shape of the on‐fault magnitude distributions and which branching faults are involved in multifault ruptures depend on how much emphasis is placed on fitting the target slip rate. In cases where the underlying data are not strong enough to warrant chasing the target slip rate, it is better to focus on the distribution of feasible results that better represents the uncertainty in the solutions imposed by the data.

We present a method to calculate landslide hazard curves along offshore margins based on size distributions of submarine landslides. The method utilizes 10 different continental margins that were mapped by high‐resolution multibeam sonar with landslide scar areas measured by a consistent Geographic Information System procedure. Statistical tests of several different probability distribution models indicate that the lognormal model is most appropriate for these siliciclastic environments, consistent with an earlier study of the U.S. Atlantic margin (Chaytor et al., 2009, https://doi.org/10.1016/j.margeo.2008.08.007). Parameter estimation is performed using the maximum likelihood technique, and confidence intervals are determined using likelihood profiles. Pairwise comparison of size distributions for the 10 margins indicates that the U.S. Atlantic and Queen Charlotte margins are different than most other margins. These margins represent end‐members, with the U.S. Atlantic margin having the highest mean scar area and the Queen Charlotte margin the lowest. We demonstrate that empirical, offshore landslide hazard curves can be developed from the landslide size distributions, if the duration of mapped landslide activity is known. This study indicates that the shape parameter of the size distribution is similar among all 10 margins, and thus, the shape of the hazard curves is also similar. Significant differences in hazard curves among the margins are therefore related to differences in mean sizes and, potentially, differences in the duration of landslide activity.

SUMMARYCombinatorial methods are used to determine the spatial distribution of earthquake magnitudes on a fault whose slip rate varies along strike. Input to the problem is a finite sample of earthquake magnitudes that span 5 kyr drawn from a truncated Pareto distribution. The primary constraints to the problem are maximum and minimum values around the target slip-rate function indicating where feasible solutions can occur. Two methods are used to determine the spatial distribution of earthquakes: integer programming and the greedy-sequential algorithm. For the integer-programming method, the binary decision vector includes all possible locations along the fault where each earthquake can occur. Once a set of solutions that satisfy the constraints is found, the cumulative slip misfit on the fault is globally minimized relative to the target slip-rate function. The greedy algorithm sequentially places earthquakes to locally optimize slip accumulation. As a case study, we calculate how earthquakes are distributed along the megathrust of the Nankai subduction zone, in which the slip rate varies significantly along strike. For both methods, the spatial distribution of magnitudes depends on slip rate, except for the largest magnitude earthquakes that span multiple sections of the fault. The greedy-sequential algorithm, previously applied to this fault (Parsons et al., 2012), tends to produce smoother spatial distributions and fewer lower magnitude earthquakes in the low slip-rate section of the fault compared to the integer-programming method. Differences in results from the two methods relate to how much emphasis is placed on minimizing the misfit to the target slip rate (integer programming) compared to finding a solution within the slip-rate constraints (greedy sequential). Specifics of the spatial distribution of magnitudes also depend on the shape of the target slip-rate function: that is, stepped at the section boundaries versus a smooth function. This study isolates the effects of slip-rate variation along a single fault in determining the spatial distribution of earthquake magnitudes, helping to better interpret results from more complex, interconnected fault systems.

An estimate of the expected earthquake rate at all possible magnitudes is needed for seismic hazard forecasts. Regional earthquake magnitude frequency distributions obey a negative exponential law (Gutenberg‐Richter), but it is unclear if individual faults do. We add three new methods to calculate long‐term California earthquake rupture rates to the existing Uniform California Earthquake Rupture Forecast version 3 efforts to assess method and parameter dependence on magnitude frequency results for individual faults. All solutions show strongly characteristic magnitude‐frequency distributions on the San Andreas and other faults, with higher rates of large earthquakes than would be expected from a Gutenberg‐Richter distribution. This is a necessary outcome that results from fitting high fault slip rates under the overall statewide earthquake rate budget. We find that input data choices can affect the nucleation magnitude‐frequency distribution shape for the San Andreas Fault; solutions are closer to a Gutenberg‐Richter distribution if the maximum magnitude allowed for earthquakes that occur away from mapped faults (background events) is raised above the consensus threshold of M = 7.6, if the moment rate for background events is reduced, or if the overall maximum magnitude is reduced from M = 8.5. We also find that participation magnitude‐frequency distribution shapes can be strongly affected by slip rate discontinuities along faults that may be artifacts related to segment boundaries.