Applying Conservation of Energy to Estimate Earthquake Frequencies from Strain Rates and Stresses

Estimating earthquake occurrence rates from the accumulation rate of seismic moment is an established tool of seismic hazard analysis. We propose an alternative, fault‐agnostic approach based on the conservation of energy: the Energy‐Conserving Seismicity Framework (ENCOS). Working in energy space has the advantage that the radiated energy is a better predictor of the damage potential of earthquake waves than the seismic moment release. In a region, ENCOS balances the stationary power available to cause earthquakes with the long‐term seismic energy release represented by the energy‐frequency distribution's first moment. Accumulation and release are connected through the average seismic efficiency, by which we mean the fraction of released energy that is converted into seismic waves. Besides measuring earthquakes in energy, ENCOS differs from moment balance essentially in that the energy accumulation rate depends on the total stress in addition to the strain rate tensor. To validate ENCOS, we exemplarily model the energy‐frequency distribution around Southern California. We estimate the energy accumulation rate due to tectonic loading assuming poroelasticity and hydrostasis. Using data from the World Stress Map and assuming the frictional limit to estimate the stress tensor, we obtain a power of 0.8 GW. The uncertainty range, 0.3–2.0 GW, originates mainly from the thickness of the seismogenic crust, the friction coefficient on preexisting faults, and models of Global Positioning System (GPS) derived strain rates. Based on a Gutenberg‐Richter magnitude‐frequency distribution, this power can be distributed over a range of energies consistent with historical earthquake rates and reasonable bounds on the seismic efficiency.


Introduction
Although the conservation of energy is a fundamental physical principle, it has not seen widespread application in seismic hazard assessment. We propose that it can be used to constrain the occurrence rates of earthquakes based on the accumulation rate of elastic energy in the crust. If the uncertainties can be sufficiently controlled, this approach would supplement the existing classes of methods now in use in the seismic hazard community: statistical inference based on seismicity catalogs and methods based on balancing the accumulated strain, inferred from geological or geodetic observations, with the strain release associated with the long-term average earthquake occurrence rates.
The first class of methods to parameterize seismicity models uses statistical inference based on seismic catalogs of the target regions (Aki, 1965;Stromeyer & Grünthal, 2015). While such methods may perform well in areas where catalogs are comprehensive, they have deficits in regions of low seismicity (Mazzotti & Gueydan, 2017) and are limited in the reliability of estimating the occurrence rates of large earthquakes.
A second class is based on the balance between accumulated seismic moment over interseismic intervals and the release of seismic moment in earthquakes. This approach incorporates long-term fault slip rates to obtain seismicity rates (e.g., Davies & Brune, 1971). From the beginning, slip rates were determined by geological observations, for example, in works by Anderson (1979), Molnar (1979), and Wesnousky (1986). In recent years, this has been complemented by the use of geodetic data (Avouac, 2015;Rollins & Avouac, 2019). Fault slip rates are useful to predict seismicity rates in regions where knowledge of the geometry and slip rates of faults is sufficiently complete. There, characterizing the seismicity of individual faults can help to create a detailed seismicity model as in case of the third California Earthquake Rupture Forecast (UCERF3) (Field et al., 2017). However, earthquakes may occur on blind faults (Talebian et al., 2004), and many shallow earthquakes do not result in a surface rupture, for example, the M W 6.8 2000 Tottori earthquake (Semmane et al., 2005) or many earthquakes of magnitude smaller than 6.6 in the Basin and Range province (dePolo, 1994). Moreover, in many areas active fault maps are still not sufficiently complete. Therefore, slip-rate based methods can often benefit from fault agnostic approaches. Kostrov (1974) recognized that the seismic moment rate in a volume could be converted to a volumetric strain rate. This method has been applied in both regional (e.g., Anderson, 1979;Bird & Liu, 2007) and global (e.g., Bird et al., 2010) applications.
We have recognized that it is possible to formulate a comprehensive new framework to link continuum mechanics to seismicity models by employing conservation of energy instead of balancing the seismic moment. The main practical difference of this new approach is that the total stress described by the stress tensor enters the stationary budget and that coupling coefficients are replaced by average seismic efficiencies.
An underlying advantage of the energy conservation approach is that the seismic energy of an event is of key interest for the assessment of ground shaking in the frequency band 1 to 10 Hz which controls most of the housing damage (Bindi et al., 2019;Kanamori et al., 2020;Picozzi et al., 2018). For instance, early estimates of the damage resulting from the M W 7.8 event in Nepal in April 2015 were overestimated as the seismic energy was considerably smaller at high frequencies because of a slow onset time (e.g., Galetzka et al., 2015). Moment-based distributions are correlated with energy-based distributions, but a perfect correlation would only occur if the ratio of seismic energy to seismic moment (or the apparent stress) was a constant (Kanamori, 1977). In reality, large variations are observed (e.g., Oth et al., 2010). In terms related directly to the seismic source, the seismic moment is determined by the amount of static deformation caused by the earthquake, while the seismic energy depends on how rapidly that deformation occurs, or the magnitude of the stresses that act to cause the slip on the fault (e.g., Deichmann, 2018).
Besides tectonic loading, the framework that we propose may be applied to all other relevant geophysical processes, such as post-glacial rebound (Nocquet et al., 2016) or erosion (Anderson, 1986), whose contribution to the seismic energy budget can be modeled. Designed with the application in cases of difficult fault identification and characterization in mind, for example, in moderate seismicity areas and areas of limited knowledge, the link is established through a fault-agnostic area source model.
We start with the description of the fundamental assumptions involved in linking the different geophysical processes to seismicity models in general. This involves energy balance, the energy-frequency relation represented by the energy-frequency density, and efficiencies. Subsequently, we describe how additional assumptions facilitate the use of interseismic velocities derived from Global Navigation Satellite System (GNSS) data, horizontal stress field orientations, faulting regimes, and elastic material parameters to quantify poroelastic geophysical processes. To illustrate the feasibility and validity, we then employ the framework and create a seismicity model for Southern California, perform a sensitivity analysis, and compare the results to observed seismicity.

Energy Balance Framework
In this section, we derive the "Energy-Conserving Seismicity Framework" (ENCOS) to link continuum mechanics to seismicity and discuss the fundamental principles and assumptions upon which it is based. ENCOS encapsulates long-term deformation processes that cause seismicity. Depending on the region of interest, different processes such as tectonic loading at plate boundary regions or weaker effects like glacial isostatic adjustment and erosion in stable continental regions (Calais et al., 2010;Steer et al., 2014;Vernant et al., 2014) may be relevant to describe seismicity. Therefore, the framework comprises a flexible method to build a seismicity model. In the following we introduce and connect the underlying principles and assumptions of the framework outlined in Figure 1. The arguably most fundamental principle of physics is energy conservation. Applied to seismicity, it implies that the combined dissipated energy of earthquakes has been converted from some sources of energy. The energy budget of individual earthquakes has been reviewed in depth by Dahlen (1977). Here, we consider the more macroscopic energy budget of the combined seismic process of many earthquakes.
Besides the fundamental energy conservation, we will make the further assumption of stationarity which is not a fundamental physical law. Assuming that earthquakes occur with steady average rates at long time scales, energy loading rates from continuum mechanics can be related to energy released in earthquakes. Consider the following stationary seismic process: The long-term continuum-mechanical process, for example, plate tectonics or glacial rebound, continuously transfers energy between different reservoirs. Among those, it continuously fills the total seismically available energy of a seismic volume of lithosphere. From this reservoir, energy is removed stepwise by the earthquake process. If the mechanical loading rate of the volume does not fluctuate significantly around a steady mean Ė, this process can be represented by the equation In Equation 1, E r is the energy radiated in an individual earthquake. The energy-frequency distribution (EFD) r(E r ) denotes the long-term average rate of earthquakes with seismically radiated energy in the infinitesimal interval [E r ,E r +dE]. The Gutenberg-Richter magnitude-frequency distribution (Gutenberg & Richter, 1944) would be an example of an EFD, if the energy radiated by the earthquake was represented by M E , the energy magnitude, that is essentially a change of variables from E r . The right hand side of the equation is the total long-term average energy radiation rate of the earthquake process (Knopoff & Kagan, 1977). The average seismic efficiency η denotes the average fraction of reservoir, that is, elastic, gravitational, and kinetic, energy released by the earthquake that is radiated as seismic waves. The remainder is dissipated essentially in friction and fault fracturing (Dahlen, 1977). We use an average seismic efficiency with the recognition that it may not be the same in every earthquake. Equation 1 is the energy conservation counterpart of similar expressions for models based on seismic moment rate conservation (Anderson, 1979;Bird & Liu, 2007).

Example: Constraints on the Gutenberg-Richter Energy-Frequency Distribution
Once the energy loading rate and the seismic efficiency are determined, Equation 1 imposes a boundary condition onto the EFD r(E r ). In principle, the shape of r(E r ) is not restricted by ENCOS; however, we will show how Equation 1 can constrain the shape of a Gutenberg-Richter relationship (GRR) (Gutenberg & Richter, 1944) employed as an EFD.
The GRR relates the long-term frequency of earthquakes with their magnitudes. Given sufficiently long observation time and large area, a region's earthquake catalog above the magnitude of completeness M c (Mignan, 2012) can generally be described by Here, N(M) is the number of observed earthquakes exceeding magnitude M in the period of observation. Equation 2 has been applied for all of the prominent magnitude scales, so we will assume that it is applicable for M E , the energy magnitude defined by Choy and Boatwright (1995). This scale is defined by the following equation: where E M0 ¼10 4.8 J is the energy of an earthquake with magnitude M E ¼0. Inserting into the GRR (2), dividing by the observation time, and deriving by energy, the power law scaling of the GRR density in energy-frequency space is obtained where β¼ 2 / 3 b. Following arguments from Knopoff and Kagan (1977), Anderson (1979), and Kagan (2002), further constraints on the EFD (4) can be identified. For a pure power law (4), the integral (1) converges at small energies only when β<1, while convergence at the infinite energy limit requires β>1. These exclusive conditions require a taper of (4) at either low or high energies. For observed natural seismicity, the b-value in moment magnitude scale lies generally below 1.2 (see, e.g., the work by Schorlemmer et al., 2005), so that β<1 and a high energy taper is required. A common choice (e.g. Main, 2000;Serra & Corral, 2017) that performs well on global data (Leonard et al., 2001) is an exponentially tapered EFD with the normalization constant C. This EFD shows power law behavior below the corner energy E c associated to the corner magnitude. Without known physical constraints to specify the taper, an exponential taper has the advantage of providing a simple and localized, that is, fast, taper while avoiding the physically unlikely discontinuous cut-off of a truncated GRR. With this parameterization of the EFD as one of a range of choices (e.g., Kagan, 2002), the integral (1) can be evaluated. The resulting algebraic equation relates C to β using Ė and η: Finally, the constant C obtained from this equation can be inserted into the EFD. The resulting EFD r γ (E r ) and its cumulative frequency distribution R γ (E r ) are functions of the four parameters Ė, η, E c , and β: Here, Γ i (−β,E r /E c ) is the upper incomplete gamma function. The role of the parameters is demonstrated in Figure 2 by varying individual parameters: The exponent β corresponds to the b-value and determines the log-log-slope. The average seismic efficiency η and the loading rate Ė control the total amount of seismicity interchangeably. Finally, the location of the taper is determined by the corner energy E c .
By relating the energy loading rate to the seismicity rates, Equation 6 imposes constraints on the EFD. Since the link is based on energy conservation and most of the seismic energy is dissipated in large earthquakes, this allows better constraining especially the recurrence rates of large earthquakes. In regions where energy loading can be inferred by geomechanical models, this method may pose a new opportunity to constrain seismicity models.

Estimating the ENCOS Parameters
The stationary energy rate conservation introduced in the previous section does not constitute an application-ready seismicity model. To implement the link described in the previous section and create an actual source model, the parameters have to be assessed. Therefore, the general steps to obtain an implementation of ENCOS can be summarized as follows: 1. Identify the key processes that control the seismically available energy loading of the target region and calculate the loading rate.
2. Choose the appropriate energy-frequency density, for example, a Gutenberg-Richter relationship (5), and relate energy loading rate to the density parameters using (1), for example, resulting in Equation 7. 3. Estimate the remaining model parameters, for example, corner energy, efficiency, and b-value, using relevant data sets of the target region.
The parameters of the EFD such as corner energy and b-value may differ between applications due to the choice of EFD. In contrast, every application of ENCOS needs an estimate of energy loading rate Ė and efficiency η . In the following sections, we consider these two parameters in detail and outline methods to estimate them.

Efficiency
Consider a geomechanical process such as tectonic loading of the crust that provides an energy reservoir driving seismicity. Let ΔE denote the energy that a single earthquake removes from the reservoir. This energy might consist of a change ΔU el in elastic and ΔU gr in gravitational potential energy and a change ΔU rot in rotational kinetic energy of the earth. While Dahlen (1977) showed that ΔE can also be expressed purely by the stress tensor and slip on the fault instead of computing the, potentially small, difference between those terms, the approach of modelling an energy reservoir in a volume considers explicitly these three contributions. The removed energy ΔE then separates into heat dissipated due to friction (E H ), fracture energy at the expanding fault tip (E G ), and energy E r radiated as seismic waves (Dahlen, 1977;Husseini, 1977): Geomechanical modeling considers the change ΔE in the energy reservoir while seismic hazard requires knowledge of the radiated energy E r . A straightforward way to transform between these two domains is the introduction of a scale factor that denotes the fraction of elastic energy that is radiated as seismic waves. Following recent convention (Bormann & Di Giacomo, 2011;Kanamori & Brodsky, 2004), we call this scale factor the seismic efficiency. Husseini (1977) separates η into the fraction η h of "available energy" after frictional dissipation and the fraction η r of the "available energy" that is finally radiated after also the fracture energy has been subtracted: Since the energies are non-negative, it follows that each efficiency takes values between 0 and 1: The radiation efficiency η r can be assessed from seismic data (Brune, 1968;Husseini, 1977;Randall, 1972;Venkataraman & Kanamori, 2004). Typical values for large earthquakes are in the range 0.25 to 1.0 (Venkataraman & Kanamori, 2004). Theoretical models (e.g., summarized by Kanamori & Brodsky, 2004) and limited observations find that η r is strongly dependent on the rupture velocity, where large rupture velocity causes high values of η r .
The efficiency η h on the other hand depends on stress tensor and the coefficient of sliding friction and cannot be inferred from seismological methods (Kanamori & Brodsky, 2004). It represents a missing link between gravitational and elastic energy ΔE and radiated energy E r . There have been efforts to directly compute the seismic efficiency by using the energy balance, for example, works by McGarr et al. (1979) and McGarr (1994). Most of these works find small values no larger than 6% for η, which indicates that η h is small. However, there are two indications that it might be considerably higher. The first is a few studies that find high values for η, for example, the largest earthquake studied by Prejean and Ellsworth (2001) with η>15% and the results of Wang (2004) with estimates of η up to 26%. The second is finding from some recent events that the dynamic friction on the fault in large earthquakes is very small (e.g., Fulton et al., 2013;Kano et al., 2006). Hence, for these large earthquakes, relatively little energy is dissipated in dynamic friction.
So far, we have considered the energy budget of single earthquakes. In the energy loading rate balance (1), the energy release of the ensemble of earthquakes has to be balanced with the total external loading. For the energy released by a finite number n of events, the energy balance reads where E r is the total radiated energy, ΔE and ΔE i are the total and individual energies drawn from the reservoir, and η i and E r,i are the efficiencies and radiated energies of the events. Hence, the average efficiency η refers to the average of the seismic efficiencies η i weighted by the energy release ΔE i and is dominated by the efficiencies of large earthquakes. Note that for the definition of the average efficiency, we used the term of an "event." This is to highlight that the energy balance includes processes like afterslip and slow slip events that are not commonly understood to be earthquakes since they do not radiate energy, but are also transient processes that release elastic energy (Avouac, 2015;Jolivet & Frank, 2020). These events are formally accounted for in the energy balance by assigning efficiency zero. In practice, this means that the average seismic efficiency can be factorized into the average seismic efficiency of earthquakes and the fraction of energy released in earthquakes: Here, η * is the average seismic efficiency of the earthquakes weighted by released energy. The fraction η eq of energy that is released in earthquakes is similar to the fraction α of seismic to total transient slip in moment balancing (Avouac, 2015). This coefficient α may on average be high (e.g., 0.8-1.0) in the upper brittle continental crust but significantly lower in deeper, hotter, or oceanic rheologies. We note, however, that while α may be an indicator for changes in η eq , the two factors are not necessarily equal: Since many transient aseismic processes are thought to occur due to stress changes, such as fluid pressure induced reductions of normal stress tensor components (Avouac, 2015;Wallace et al., 2017), and the change of elastic energy is sensitive to the stress tensor as shown, for example, in the following section, it is well possible that the impact of aseismic transients on the energy budget differs from its impact on the moment balance.
Due to neither the total earthquake energy nor the radiated energy being easily measured and since energy balance has not been in focus of research in the last decades, empirical data on earthquake efficiency are rare.
In particular, no comprehensive study of the efficiency for different settings of seismicity exists. Therefore, the efficiency remains one of the key uncertainties within ENCOS. This will be addressed in section 5.

Example: Poroelastic Energy in a Depth-Symmetric Critically Stressed Crust in Hydrostatic Equilibrium
The energy potentially available to produce earthquakes is the surplus of reversible deformation, rotational kinetic, and gravitational energy compared to the maximally relaxed state. In this section, we exemplarily consider regions where changes in gravitational and rotational energy are negligible, which we assume might be the case in tectonic settings dominated by strike-slip earthquakes. In such regions, the energy available from within a volume V to produce earthquakes is determined by integrating the reversible deformation energy density w over V. Consequently, a stationary energy loading rate is obtained by integrating the time derivative _ w over the volume: To demonstrate the application of ENCOS we choose linear poroelasticity as a simple geomechanical model. In Appendix A we show how, following Biot (1941Biot ( , 1962 and Wang (2000), the poroelastic power density can be written in terms of the strain rate tensor _ ϵ, the stress tensor σ, and the pore pressure p and its derivative using the Biot-Willis coefficient α and the constrained specific storage coefficient 1/M: This formulation of the energy loading rate density highlights a difference between the ENCOS and seismic moment balancing approaches: While the latter depend only on the strain rate tensor (e.g., Avouac, 2015;Bird et al., 2010), the elastic energy balance depends additionally on the static stress tensor or equivalently the static strain tensor. One consequence is that the relative orientation of the strain rate and stress tensor is important for the energy balance. In Appendix B, we show that in a two-dimensional system, the trace of the product can be written as tr σ _ ϵ ð Þ ¼ S Hmax ðcos 2 ðθÞϵ 11 þ sin 2 ðθÞϵ 22 ÞþS hmin ðsin 2 ðθÞϵ 11 þ cos 2 ðθÞϵ 22 Þ where θ is the angle between the principal coordinate systems of the two tensors and S Hmax and S hmin are the maximum and minimum horizontal stress, respectively. To highlight the impact, in the hypothetical case that the angle between these systems were 45°and p¼0 and we would observe pure shear (ϵ 11 ¼−ϵ 22 ), the elastic power would be zero while deformation would be visible.
Equation 16 allows for a range of helpful approximations that we will consider in the remainder of this section. We start with the assumption that stress and strain rate tensors are depth translation invariant except for the three principal stress magnitudes of the stress tensor. This assumption is made for strike-slip systems of depth-symmetric fault geometries which are driven from the side, the latter being an approach taken also in existing models (e.g., Bird, 2009). Then, the strain rate tensor is determined by the surface velocity field v ! by where we have chosen a Euclidean coordinate system and strain to be positive for extension. The surface velocity field can in turn be estimated from geodetic observation techniques.
For the stress tensor, data on the orientation of S Hmax are provided in the World Stress Map (WSM) database while magnitudes are seldomly available (Heidbach et al., 2018). A simple analytical approximation of the 10.1029/2020JB020186

Journal of Geophysical Research: Solid Earth
stress magnitudes is the concept of the critically stressed crust (Jaeger et al., 2007;Sibson, 1974). We employ it by assuming that the local stress tensor is determined by nearby, possibly unknown, faults characterized by cohesionless Mohr-Coulomb fracture criterion with friction coefficient μ. The stress tensor is then determined by the gravitational load, the local style of faulting, the friction coefficient μ, and the S Hmax orientation. In Appendix C, we describe in detail how to compute the stress tensor.
A final approximation we employ is to assume hydrostasis, that is, _ p ¼ 0. This is motivated by the observation that while pore pressure variations play an important role in the dynamics of faulting (e.g., Avouac, 2015), pore pressure in the bulk material, which dominates the balance (15), is often found to be hydrostatic (Townend & Zoback, 2000).
In the depth-invariant approximation, the depth dependence of the energy density _ wð x ! Þ separates from the horizontal components and is uniquely captured by the depth of the volume storing elastic energy for earthquakes, that is, the variable thickness of the seismogenic crust z max . Hence, the depth dependence of _ w can be collapsed into an area density _ w a ðx; yÞ of energy loading rate that will depend quadratically on the depth of the contributing volume: The exact manifestation of this depth depends on the geological setting and the geometry. Our working hypothesis is that z max is related to the brittle-ductile transition zone which may decouple the brittle crust from greater depths and release strain plastically. From the seismological perspective, this depth may be estimated from the depth distribution of earthquakes in the region.

Testing ENCOS in Southern California
With the complete description of ENCOS at hand, its essential validity remains to be tested. To demonstrate its application, we select Southern California surrounding the San Andreas fault, a region characterized by tectonic loading (e.g., Anderson, 1971) and well-monitored seismicity. These properties qualify the region as a simple feasibility test case.
We choose our region of interest (ROI) with the depth-symmetric two-dimensional poroelastic approach and a good data coverage in mind. It spans from the end of the creeping section to the Gulf of California and from the Walker Lane into the Pacific Ocean. In an oblique Mercator projection (Snyder, 1987)-minimizing distortion across the ROI-a rectangular coordinate raster is chosen. The resulting ROI is shown in Figure 3. To avoid aliasing, we use an 800 m resolution in this raster to match the maximum resolution attainable in any region of any of the data sets used. Before combining the fields obtained from different data sets, we smooth each with a second-order Butterworth low-pass filter of data-driven bandwidth λ s . Before smoothing, we apply a two-dimensional cosine taper configured such that the study area data remain unchanged. The outer box in Figure 3 shows this margin for the present example. After smoothing, we renormalize the fields to the root mean square of the unprocessed fields to prevent an energy reduction of the individual fields purely due to the smoothing process.
The smoothing bandwidth λ s is an estimate of the maximum resolution achievable from the data. It is determined by the largest resolution-inhibiting effect we estimated for all spatial data sets used. For our test area in Southern California, this is the spatial distribution of the World Stress Map (WSM) data (Heidbach et al., 2016). We estimate the resolution by viewing the data points as spatially distributed sensors. Based on the concept of the seismic array response function (Capon, 1969;Rost & Thomas, 2002), we estimate the spatial bandwidth in which the spatial configuration can detect patterns. We obtain a band from 32 to 240 km. From the band, we use twice the logarithmic mean of the bounds, λ s ¼130 km, as preferred value and consider the latter in the uncertainty analysis.

Stress and Strain
For the estimation of the S Hmax orientation, data records are taken from the WSM database release 2016 (Heidbach et al., 2016(Heidbach et al., , 2018) based on information from earthquake focal mechanisms, borehole and geological data. We only use data records with A-C quality except those earthquake focal mechanisms that are labeled as possible plate boundary events according to Heidbach et al. (2010) or that are older than 1976. Furthermore, we excluded three data records from hydraulic fracturing tests and borehole breakouts, respectively, where detailed information for a reproducible quality assignment is missing.
To estimate the local faulting type, we use Nadaraya-Watson interpolation (Nadaraya, 1964;Watson, 1964) with a Gaussian kernel parameterized by bandwidth λ g . The interpolated faulting type is based on the WSM data records that also provide the information on the stress regime. In areas with ambiguous faulting type, the stress magnitudes for the intermixed faulting regimes are linearly interpolated. For the frictional limit derived stress magnitude for strike-slip regimes (C6), we use κ¼0.5 (see Appendix C). For all faulting types, rock density is estimated as homogeneous bulk property from the CRUST1.0 model (Laske et al., 2013) inside the study area.  Felzer and Cao (2008) with M W ≥ 6.5. Circle size represents energy magnitude M e (Bormann & Di Giacomo, 2011) for the IRIS catalog, else M W . In Figure 9, the earthquakes within the inner area are compared to the ENCOS seismicity model of the inner area. We include in that comparison also the 1857 Fort Tejon earthquake, whose epicenter lies just outside the inner area, since its rupture spans into the inner area. Topography is adapted from ETOPO1 model (Amante & Eakins, 2009), coastlines are taken from the GSHHS database (Wessel & Smith, 1996), and fault traces from the UCERF3 dataset (Field et al., 2013).
Finally, the static friction coefficient μ needs to be estimated. For friction coefficients in California, a variety of different results have been published (Brune & Thatcher, 2002;Scholz, 2006). Values inferred from individual borehole measurements range from μ¼0.1 in the gouge zone of the southern part of the San Andreas creeping section (Carpenter et al., 2015) to larger values 0.6 ≤ μ ≤ 1.0 in boreholes surrounding the San Andreas fault (Zoback & Healy, 1992). Due to the heterogeneity of the borehole results and the sparse measurement coverage of the ROI, we use a value μ¼0.4 between the measurement extremes and later consider the extremes in the uncertainty analysis (see Appendix E). The resulting stress tensor field is shown in Figure 4. Since due to Equations C1 and C3 the stress tensor is linear in z, we present its depth derivative.

Strain Rate Tensor
For simplicity and due to a lack of deformation-indicating data at depth, we use the surface approximation (18). We compute these surface strain rates directly from interseismic GPS velocities evaluated using the MIDAS method (Blewitt et al., 2016(Blewitt et al., , 2018 for our central estimate and for comparison in the uncertainty analysis using the global strain rate model (GSRM) version 2 (Kreemer et al., 2014). In Appendix D, we describe the preprocessing of the MIDAS velocity data to remove outliers and stations whose results are dominated by anthropogenic signals that mask the tectonic signal relevant to seismicity. The local computation is performed by nearest neighbor interpolation of the GPS velocities followed by spatial differentiation for use in Equation 18 and subsequent smoothing of the tensor field. The central estimate is displayed in Figure 5.

Total Power and its Sensitivity
Again using the Nadaraya-Watson interpolator with Gaussian kernel and bandwidth λ g , we estimate the local depth z max of the loaded volume from the maximum fault depth given by the UCERF3 Fault Model FM3.1 (Field et al., 2013(Field et al., , 2015. This goes along with observations that for Southern California, the brittle-ductile transition zone estimated from heat flow relates to the depth of the seismogenic zone (Bonner et al., 2003;Hauksson, 2011). The depth field obtained is shown in Figure 6. Following Equation 19, the depth-integrated energy loading rate density _ w a is calculated and shown in Figure 7. The total poroelastic energy loading rate then follows by integrating _ w a over the study area: The value of 0.8 GW, in the order of magnitude of a large electrical power plant, is calculated for the preferred scenario with one significant digit. However, as mentioned in the previous section, the parameters, and some of the methods, are subject to considerable uncertainty.
To quantify the results' uncertainties and the impact of the respective uncertainties in the input data, forward uncertainty propagation is applied in two steps. First, all physical constants involved in the calculations and all parameters of the algorithms used are varied within reasonable bounds. The specific values of these bounds are discussed in Appendix E. Where impact could be considered insignificant a priori, constants were not varied. In a second step, a Monte Carlo approach was used to vary the input data of the WSM according to the uncertainties of the individual data points. This analysis was done using multivariate  Figure 8 and yields a range of for the elastic energy loading rate. The most obvious finding is that the friction coefficient μ and the maximum depth z max dominate the estimate of the uncertainty. Granted that we analyzed a rare extreme case of homogeneous z max set to extreme values obtained from Nazareth and Hauksson (2004) and that _ w a in Equation 19 scales quadratically with z max , the high impact of z max is expected for our method.
The friction coefficient μ is a near-linear constituent of the stress tensor-and consequently the energy loading rate-due to the failure criterion. Since we cannot safely discriminate between fracture of intact rock, failure of preexisting faults, or the precise mixture of both when it comes to the prevalent rupture scenarios in the volume, the friction coefficient has a high plausibility range and thereby is a key driver of the uncertainty.
The smoothing bandwidth λ s has a small negative impact. This may indicate its correctness: For larger λ s , the energy inherent to high modes of the fields is redistributed to lower modes. Therefore, not much is expected to change. For smaller λ s , that is, higher detail, it seems possible that the densely distributed GPS stations measure small-scale perturbations of the large-scale strain rate field. Due to the lower density of stress tensor proxies, any stress tensor structure corresponding to these strain rate details cannot be detected. Hence, on  (Field et al., 2013). They are used to estimate the interpolated volume depth. Near-surface faults shallower than 8 km and faults close to deeper reaching faults (distance smaller than depth of the deep reaching fault) are both excluded from the interpolation.

Journal of Geophysical Research: Solid Earth
small scales the strain rate tensor could deviate from an energy-maximizing orientation relative to the stress tensor. This would create patches of smaller energy loading rate compared to the large-scale estimate. A sufficiently large smoothing bandwidth would circumvent such an effect. The exact nature of this mechanism might be interesting for future study. For the purpose of this work, the impact of the smoothing bandwidth is a second-order uncertainty.

Journal of Geophysical Research: Solid Earth
The remaining quantities are rather insignificant to the total uncertainty of the energy loading rate. In particular this includes the Biot-Willis coefficient α, the only poroelastic parameter remaining in Equation 16 for hydrostatic conditions.

Comparing to Observed Seismicity
Now that the seismic energy loading rate of the source volume has been calculated, the result can be compared to earthquake catalogs to investigate whether the model with its linear elastic energy loading rate budget is consistent with the observations. Besides estimates of the model parameters, this comparison requires catalogs of earthquake energies. In the ROI, two types of catalogs are available: historical and instrumental catalogs. A comprehensive catalog of historical earthquakes has been compiled by Felzer and Cao (2008) for the UCERF2 model. Additionally, starting 1990 the area is covered by the global energy catalog of the Incorporated Research Institutions for Seismology (IRIS) (Convers & Newman, 2011;IRIS DMC, 2013), and Yang et al. (2012) have computed a refined focal mechanism catalog spanning the years 1981 to 2019. Only the IRIS catalog measures earthquakes in energy, and it does so only for earthquakes of magnitude M W ≥ 6.0. This is due to the sensitivity of the estimation of the radiated energy to directivity and attenuation effects which make E r difficult to determine (Deichmann, 2018;Kanamori et al., 2020;von Specht et al., 2019). Recent efforts to routinely calculate energy magnitudes might increase earthquake energy data availability in the future (Strollo et al., 2020). The other two catalogs are given in moment magnitude, which is largely unaffected by attenuation and directivity effects (Deichmann, 2018), and is readily available through stable methods (Bormann & Di Giacomo, 2011). As noted in section 1, there is no perfect relation between the moment magnitude and the radiated energy since the former is a kinematic and energy a dynamic measure. However, the empirical correlations between the seismic moment and the radiated energy can be used to convert between the two on average. For the introduction of the moment magnitude, Kanamori (1977) used a relation first introduced by Gutenberg (1956): For individual earthquakes, the radiated energies scatter significantly about this empirical equation (Bormann & Di Giacomo, 2011;Picozzi et al., 2018). On average, however, relation (22) is supported by both recent empirical (Kanamori et al., 1993(Kanamori et al., , 2020Picozzi et al., 2018) and theoretical studies (Deichmann, 2018) for earthquakes larger than M W ¼4.0. Therefore, we include in our comparison the moment magnitude catalogs by Felzer and Cao (2008) and Yang et al. (2012) in addition to the earthquake energy catalog by IRIS DMC (2013) and compute earthquake energy from M W using Equation 22. We retain only earthquakes that have been located within the ROI as shown in Figure 3, except for the 1857 Fort Tejon earthquake whose epicenter lies just outside the ROI but whose rupture extends significantly into the ROI. The two moment magnitude catalogs overlap partly from 1981 onward. To prevent double counting of those earthquakes, we remove the overlapping part of the historical catalog. The two instrumental catalogs overlap over the full time span of the IRIS catalog. Nevertheless, we retain all earthquakes of both catalogs to illustrate the differences in energy estimates of both methods, direct and using Equation 22, and to assert that the empirical relation does not introduce a significant bias in our comparison. Following the assessment of Felzer and Cao (2008), we assume a magnitude of completeness of M W 6.5 for the historical catalog from 1850 onward which retains 22 used earthquakes. Figure 9 shows the high-energy range of the empirical cumulative distribution , which has also been used to convert the moment magnitude catalogs to earthquake energies. function-the cumulative number of earthquakes divided by the catalog time span-for all three catalogs. Since the uncertainty of the magnitude estimates of the historical catalog is diverse and mostly significantly large, we show also the uncertainties of the magnitude estimates for this catalog. The instrumental catalogs are complete in the shown magnitude range (Schorlemmer & Woessner, 2008) with 10 and 5 earthquakes, respectively.
In addition to the catalogs, Figure 9 contains three different parameterizations of the EFD (7) for comparison. The parameters of the EFDs are chosen to reflect the range of plausible values. For the average seismic efficiency η, we use the value range 0.2% to 2.8% obtained by McGarr (1994) for mine collapses using the driving gravitational loading. Additionally, we show one parameterization with a higher η of 6.0% following the arguments presented in section 3.1. For β, we select a range of values according to b values from b¼1 down to b¼0.7, compatible with the range of results by Schorlemmer et al. (2005).
Finally, the corner energy E c remains to be selected. In our ROI we assume that the largest earthquakes will be on the San Andreas fault system. The northern edge of the region is at the southern end of the creeping section of the San Andreas fault, so we assume that is the northern limit of rupture of the maximum earthquake. The southern limit is unclear: A rupture starting in the north could follow either the San Andreas or the San Jacinto fault to nearly the southern end of our region. From there, either of these ruptures would need to jump over a significant step to join the Imperial-Cerro Prieto fault to the southern limit of the ROI. Nonetheless, we follow recent suggestions implying that such a jump could occur (Dorsett et al., 2019) and thus measure the rupture over the full length. This gives even larger ruptures spanning beyond the ROI. For the fault section contained within the map with a length of roughly 700 km, we apply the scaling model M3 of Anderson et al. (2017) both with slip rate between 20 and 50 mm/yr and without slip rate. This corresponds to median M W between 7.8 and 8.0 and an uncertainty of 0.2 magnitude units. We note that this is larger than any earthquakes of the past ∼600 years (Scharer et al., 2014(Scharer et al., , 2017, but such ruptures are allowed by UCERF3, so we consider this a reasonable maximum energy. Following the arguments in Appendix F, we choose our central estimate of E c to be 0.2 magnitude units below the largest rupture scenario from the scaling model. Considering the uncertainty of the scaling relationship, we choose M W between 7.8 and 8.0 as estimates for E c .

Discussion
From Figure 9 it follows that within the limits imposed by the uncertainties of the individual parameters, the cumulative EFDs parameterized using ENCOS reasonably enclose the empirical cEFDs of the catalogs. Hence, within those limits, we have found ENCOS to be fundamentally plausible. While we have been able to reproduce observed seismicity using only geodetic and geomechanical data in combination with rupture size relations and a four-parameterical model, the uncertainties are still large. To constrain the uncertainties of the efficiency, it is vital to establish the dominant faulting dynamics. In cases where the enlargement of the damage zone contributes significantly to the energy budget, the seismic efficiency may be low in the order of 1%, while negligible new damage and low friction on weak faults could potentially lead to larger seismic efficiencies.
In addition to the uncertainty whose effects are shown in Figure 9, nearly one order of magnitude of energy loading rate uncertainty remains (see Figure 8). This uncertainty is mainly controlled by two parameters: the friction coefficient μ and the depth z max of the bulk volume. Observe that the static friction coefficient used to estimate the stress magnitude from the frictional limit might be correlated with the dynamic friction coefficient that controls the seismic efficiency. Hence, any increase in knowledge about the frictional process that concerns both static and dynamic friction may improve the estimate of both the elastic loading rate and the seismic efficiency. In other words investigating the friction coefficient may benefit the assessment of the seismic efficiency and vice versa.
Finiteness of the energy rate balance requires the existence of a sufficient tail decay of the EFD, for example, some sort of upper-bound cut-off, when applying ENCOS. For the Gutenberg-Richter EFD we chose an exponential tail, which introduces only one additional parameter, the corner energy E c , and provides a smooth yet rapid tapering. While this matches our expectation that earthquake energies are not discontinuously capped at the largest magnitude deemed plausible, this choice is not based on mechanical first principles. Hence, the high magnitude shape presents a large uncertainty of the model.
In principal, ENCOS is fault agnostic. In our exemplary comparison, we have used knowledge about the fault system in one occasion, to estimate the size of the largest plausible earthquake. This specific way of parameterizing an EFD is not a fundamental requirement of ENCOS; the framework imposes constraints on any EFD supplied.
In our particular application to a region dominated by transform faults, the focus on the elastic energy with the assumptions of the frictional limit and the two-dimensional approximation of strain rate and stress tensor fields seem to work well for the bulk volume. However, there are a range of unknowns, and some adjustments required for other areas that are beyond the scope of this work. For different regions, such as subduction zones, the two-dimensional approximation may fail, estimation from GNSS data may be more difficult, gravitational potential and rotational energy would have to be included, or the frictional limit may not be suitable. There, more elaborate methods such as 3-D and 4-D modeling of the strain rate and stress field Reiter & Heidbach, 2014) and additional data sets would have to be applied. Generally, the use of long-term dynamic modelling and the inclusion of data indicating long-term deformation may also be helpful in assessing the stationarity of the deformation rate, an open question ENCOS shares with moment balancing (Avouac, 2015;Holschneider et al., 2014). In some cases, nonlinear rheology may be important for the geomechanical analysis (Peña et al., 2019), although in case the nonlinear rheologies are irreversible, the elastic approximation can present an upper bound on the available power. Furthermore, in some regions, aseismic events may play a significant role in the release of seismic moment (Avouac, 2015;Hirose et al., 1999;Jolivet & Frank, 2020;Wallace et al., 2017) without radiating a significant amount of energy, which may be another significant contribution to the average efficiency. An interesting issue arising is how the seismic moment coupling coefficient compares to the fraction of aseismic strain energy release, the latter depending additionally on the stress tensor.
To apply ENCOS, key processes controlling the loading of the seismic energy reservoir need to be identified. In our test case around Southern California, this reduces to the strain energy loading due to tectonic loading. In other regions such as low-strain areas, elastic tectonic loading may not be the single dominant contribution to the energy reservoir. There, other processes such as crustal flexure due to glacial isostatic adjustment (Nocquet et al., 2016) or redistribution of vertical crustal loading due to erosion (Anderson, 1986) may be relevant contributions to the stationary energy budget. Furthermore, we expect that ENCOS can be used to test the hypothesis of Calais et al. (2016) which says that in stable continental regions transient changes of the rock strength rather than the tectonic loading rate control the seismicity rate.
In this work, we have used earthquake catalogs only for validation purposes. In application scenarios, such as the formulation of probabilistic seismic hazard assessment (Cornell, 1968;Woessner et al., 2015) in energy space, ENCOS may be combined with catalog data by means of statistical methods to estimate EFD parameters. This utilizes all available data and allows the self-consistent quantification of statistical uncertainties.
Finally, we emphasize both similarities and differences between ENCOS and the approaches based on balancing seismic moment. With our application to Southern California, we have chosen an area where seismicity is described successfully by a range of detailed models employing moment conservation and we have made use of the correlation between radiated energy and seismic moment. Hence, this example shows that to first order ENCOS and moment-based methods can be used equivalently in these tectonic settings. A technical aspect of this is that energy opens a new viewpoint on stationary seismicity and can be a potentially useful modelling target for the aforementioned nontectonic processes. Besides this, the use of energy conservation has two main fundamental differences to moment conservation: First, on the available side of the budget, ENCOS takes into account the stress tensor while moment-based methods use only the strain rate tensor (e.g., Bird et al., 2010). Coupling coefficients can calibrate both methods if the product of strain rate and stress tensor is homogeneous. In such homogeneous cases, an independent inversion of strain rate and stress tensors might introduce artificial uncertainties, and also coupled inversions of both tensors taking into consideration the physical dependencies in ENCOS might tend to yield larger uncertainties in estimating the budget than a moment balancing approach. However, by virtue of Equation 19, there can be differences when the orientation between strain rate and stress tensor is indeed heterogeneous across the region of 10.1029/2020JB020186

Journal of Geophysical Research: Solid Earth
interest. Therefore, when seismicity models are formulated in seismic moment or Equation 22 is used to convert from moment catalogs, ENCOS can at least provide an additional method of validation due to its sensitivity to these inhomogeneities.
Second, on the release side of the budget ENCOS measures earthquakes in energy which has been argued to capture the shaking potential better than seismic moment (Bindi et al., 2019;Bormann & Di Giacomo, 2011;von Specht et al., 2019). While recent analyses have confirmed that the correlation between seismic moment and radiated energy is fairly stable on average in Japan (Kanamori et al., 2020), fluctuations remain (Oth et al., 2010), and measuring earthquakes directly in energy nevertheless removes this layer of fluctuations from the hazard estimation.

Conclusions and Outlook
We presented the Energy-Conserving Seismicity Framework (ENCOS), a fault-agnostic method to combine strain rate and stress data into a stationary regional seismicity model by means of energy conservation. Employing an elementary energy-frequency distribution (EFD) and paleoseismic estimates of the maximum earthquake energy, the model estimates the rates of large earthquakes in Southern California. Our estimate of the energy loading rate is 0.8 GW, and our sensitivity tests suggest a range from 0.3 to 2.0 GW. Within present uncertainties of the EFD parameters, the results are compatible with empirical long-term rates of historical and instrumental catalogs.
ENCOS transfers the structure of fault-agnostic stationary moment rate balance (MRB) to energy space. Since both energy and seismic moment are predominately released in large earthquakes, ENCOS presents a method to estimate the rates of large earthquakes similarly to fault-agnostic methods based on conservation of moment. Due to the highly correlated nature of radiated energy and seismic moment, ENCOS and MRB are interchangeable to first order in many use cases. This can potentially present an opportunity in cases where energy as an interdisciplinary language is more accessible to a problem description than seismic moment.
The fundamental distinction between ENCOS and MRB is that the energy rate balance has an additional dependence on the stress tensor. Potential differences in moment and energy production seem possible if the relative orientation of strain rate and stress tensors is inhomogeneous, and differences in moment and energy release are at least locally possible due to variations in coseismic or aseismic stress changes. Noting the recent arguments that energy would be a more suitable scalar representation of earthquake hazard than moment magnitude since radiated energy is a better estimator of ground motion (Bormann & Di Giacomo, 2011;Bindi et al., 2019), these subtle differences make conservation of energy based seismicity models both an interesting tool and subject for further study.

Appendix A: Poroelastic Power Density
The linear poroelastic energy density w can be derived from the expressions given by Biot (1962): Here, ζ is the increment in water content, 1/H the compressibility, and I ϵ 1 and I ϵ 2 the invariants of the strain tensor ϵ defined as Using the relations given in the same work, H¼λ c +2G, λ c ¼λ+α 2 M, C¼αM, and p¼−αMε+Mζ, the energy density resolves to 10.1029/2020JB020186

Journal of Geophysical Research: Solid Earth
Here p is the pore pressure, M the inverse constrained specific storage coefficient, G and λ the Lamé constants, and ε and I ϵ 2 the first and second invariants of the strain tensor. G is also called the shear modulus. By time differentiation, the corresponding energy loading rate density, or power density, follows: The term _ ε is the first invariant of the strain rate tensor _ ϵ. The time derivative of the second invariant contains both the strain tensor and the strain rate tensor: The constitutive equation for isotropic media relates stress tensor σ and strain tensor ϵ with Poisson's ratio ν, the Biot-Willis coefficient α, and the drained bulk modulus K (Biot, 1941;Wang, 2000). In the convention we use following Wang (2000), the stresses are negative for compression, and the volumetric strain is positive in expansion, that is, defined by the following equation with the transformation u ! of relaxed to strained material space: Using Equation A7 and relations between the poroelastic constants (Kümpel, 1991) λ can be rewritten (see supporting information Text S1) into the form (16) which depends on the strain rate and stress tensor and is independent of purely elastic parameters.

Appendix B: Power Density: Strain Rate-Stress Product in 2-D
The poroelastic power density (16) depends on the trace of the product of the strain rate tensor _ ϵ and the stress tensor σ. In general, this product depends on the relative orientation of the principal coordinates of the two tensors. We illustrate this with a simple example of two-dimensional strain rate and stress tensors. The trace of the product of both tensors can be computed in the stress tensor's principal coordinate system by rotating the strain rate tensor from its principal coordinate system where θ is the angle between both coordinate systems. This equation evaluates to (17).

Appendix C: Stress Tensor in Frictional Limit
Here we describe in detail the estimation of the stress tensor in the frictional limit. First, we assume that gravitational vertical loading with rock density ρ r and gravitational acceleration g defines one of the principal stresses (as, e.g., in Reiter & Heidbach, 2014;Zang & Stephansson, 2010). It follows that the other two principal stresses are the maximum and minimum absolute horizontal stresses S Hmax and S hmin , respectively. Then, in a local coordinate system rotated in the xy plane such that the main principal horizontal stress axis coincides with the x axis, the stress tensor is of the form where z increases with depth. The tensor can then be rotated back to the initial system.
Having obtained the orientation of the principal stresses and the magnitude of the vertical principal stress, the magnitudes of the horizontal principal stresses S Hmax and S hmin remain to be determined using the concept of the critically stressed crust (Jaeger et al., 2007;Sibson, 1974). We start with including the drained pore pressure corresponding to a porously connected lithosphere. This reduces the normal stress that acts on the faults and controls the strength of the shear-resisting static friction, while the shear stress is not affected (King Hubbert & Rubey, 1959). Thus, we have to use the effective vertical stress to derive the failure condition (Shebalin & Narteau, 2017;Sibson, 1974 Assuming that the vertical stress magnitude lies within the range of horizontal stress with 0 ≤ κ ≤ 1, applying the Mohr-Coulomb failure criterion leads to the following difference between the two horizontal principal stresses (Sibson, 1974): Here R′ depends on the friction coefficient μ (Sibson, 1974;Shebalin & Narteau, 2017) and the parameter κ scales between the three faulting regimes (normal, strike-slip, and thrust faulting). In the thrust faulting (TF) regime, the vertical stress is the least principal stress (jS TF Hmax j > jS TF hmin j ≥ jσ v j). We approximate S TF hmin with σ v (κ¼0): In the normal faulting (NF) regime, the vertical stress is the largest principal stress (jσ v j ≥ jS NF Hmax j > jS NF hmin j). We approximate S NF Hmax with σ v (κ¼1): In the strike-slip faulting (SSF) regime, the vertical stress is intermediate principal stress (jS SSF Hmax j > jσ v j > j S SSF hmin j):

10.1029/2020JB020186
Journal of Geophysical Research: Solid Earth Here, 0<κ<1 scales between the end members for thrust and normal faulting.

Appendix D: Southern California: Outliers in GPS Velocities
Since its large number of stations and the automatic processing strategy fits well with ENCOS, we use the GPS velocity solutions from the dataset by Blewitt et al. (2018) obtained using the MIDAS algorithm (Blewitt et al., 2016). Its automatic and robust classification make it a promising candidate exactly for regions where dense GPS networks can deliver precise information on the strain rate fields. However, even MIDAS requires a minimum time series quality, such as a maximum fraction of jumps. If a time series does not adhere to this data quality, MIDAS may fail to reproduce the long-term velocity relevant to straining of the crust. Since we interpolate the point measurements of the GPS time series, single such outlier stations in sparsely covered regions may impact the strain rate field on a large area. Therefore, these stations may significantly impact the estimated energy loading rate.
To reduce the possible impact of outlier results, we test the impact each station has on the strain rate field and test the GPS stations that have the highest impact. First, we remove stations PLO1, ORE2, ORE3, and COMA that have obvious large displacement artifacts due to non-tectonic origins. Additionally, we remove FCTF in Los Angeles which is dominated by uplift. For all other nodes whose individual removal changes the RMS of the strain rate field by more than 1.0%, we evaluate the time series of the station for possible artifacts that could justify classifying the station as an outlier. We then remove the following stations' MIDAS GPS velocities from the data set: • Station WLHG has a single station impact of close to 16% on the strain rate field RMS. At the same time, the velocity is anomalous by 4 mm/yr in east direction compared to surrounding stations. This matches over the observation period of 4 years by a discontinuous step of 4 to 6 cm with missing data in the first of half of 2017 not considered in the MIDAS trend. • Station NVLX with an impact of 7% shows dominant signals of subsidence caused by anthropogenic depletion of aquifers in the Mexicali valley (Glowacka et al., 2010). • Station BFSH with an impact of 4% also shows a nonlinear time series with nonperiodic undulations, fits the MIDAS trend poorly, and differs in velocity strongly from station CAIN just 2.1 km away and showing a very linear trend. • Station CAFP and CAHA with an impact of 2% and LEMA and MULN with an impact of 1% lie within the water-depleting San Joaquin-Huron area (Sneed et al., 2018) and show dominant subsidence behavior possibly masking horizontal tectonic strain rates. • Station BUEG with an impact of 1% contains an undetected step.
The resulting GPS velocity data set is not guaranteed to be cleaned of outliers, but individual outlier impact is bound by 1% on strain rate field RMS. This means that data points either align well with their local neighboring stations, that is, are consistent under the assumption of a smooth strain rate field, or that they lie in densely sampled regions where their impact is limited.

Appendix E: Southern California: Estimated Uncertainties of Model Parameters
For the tornado plot Figure 8, we varied the parameters of the geophysical model and measured the response of the energy loading rate. For the variation, we estimated uncertainty intervals within which the individual parameters may reside for Southern California. This is a rather crude estimate but should give a first estimate of the impact of the individual parameters and the order of the uncertainty. We characterize the parameter bounds as follows: • For the model volume depth z max ð x ! Þ, we choose global maximum and minimum bounds for all model cells. For the maximum depth, we choose 20 km corresponding to the maximum seismogenic depth in the area obtained by Nazareth and Hauksson (2004). For the minimum z max , we choose z max ¼11 km close to the minimum z max of the preferred model.

Journal of Geophysical Research: Solid Earth
• For the purpose of computing the elastic energy loading rate using the frictional limit, it is important that the friction coefficient reproduces the volume-average stress magnitudes. In the bulk of intraplate regions, observations agree well with a critical stress state with μ>0.6 (Townend & Zoback, 2000). However, the San Andreas fault is a geometrically significant geological feature of the study area and geomechanical models have shown that large weak faults can influence the studied stress field Hergert et al., 2015). While the low friction coefficient of the creeping section (Carpenter et al., 2015) may not be representative of the whole San Andreas fault (Scholz, 2006), a weak San Andreas fault nevertheless represents an important possible extreme case of the state of stress in Southern California. Hence, to explore the range of the friction coefficient μ, we choose a maximum value of 0.85 corresponding to intermediate pressure direct contact sliding (Byerlee, 1978), and a minimum value of 0.08 corresponding to the low results of Fulton et al. (2013) and the values measured from the SAFOD borehole (Carpenter et al., 2015). • We obtained ρ r bounds from the standard deviation of the CRUST1.0 model over the study area (Laske et al., 2013). • For the smoothing bandwidth λ s , we considered the interval of length scales in which the combined data sets are sensitive, that is, a smoothing bandwidth from 32 to 240 km. Sampling this range allows us to quantify the impact of our ad hoc choice of bandwidth within this range and of yet unspecified possible inhomogeneities in the combination of data sets. • We vary the grid bandwidth λ g within a factor three from 1.2 to 3.2 km to assure that it has no significant impact. • We vary the number of data points N wsm required in each search of the WSM interpolation (Ziegler & Heidbach, 2017) between 5 and 10. Additionally, the kernel is varied between Gaussian, linear decreasing, and uniform distance weighting. • Due to a lack of data, we allow maximum uncertainty for the Biot coefficient and test the whole value range from 0 to 1 with a default of 0.5. • Strain rate tensors are estimated both directly from GPS data as well as interpolated from the GSRM (Kreemer et al., 2014).
Given these parameter bound estimates, we then vary individual parameters within these bounds and capture the extreme variations in the tornado plot Figure 8.

Appendix F: Loading Rate Exceeding Multiples of E c
Here, we consider the role of the corner energy E c of the Gutenberg-Richter type model (7). For energies exceeding the corner energy, the occurrence rates of earthquakes rapidly converges to zero. Due to the smoothness of the convergence, multiple ways exist to determine the corner energy from an upper plausible limit E m , for example, obtained from geology. One possible way is to require the occurrence rate of earthquakes exceeding this E m to be less than a rate R m . This rate could be determined such that the occurrence probably of earthquakes exceeding E m is negligible over a design time.
Another way which we employ in this work is to consider the combined power P ex of earthquakes exceeding energy E m . One can require this power to be insignificant compared to the total elastic loading rate Ė. If the corner energy is chosen to be the relative power Π ex exceeding E m amounts to We show this excess power in Figure F1. For typical values of b and a multiple α¼4, the model places a fraction of about 1% of the energy to be released in earthquakes exceeding E c . In other words, for the energy leaking into earthquakes exceeding implausibly large earthquakes to be less than 1%, one can choose E c ¼0.25E m .
On magnitude scales, this multiple α corresponds to choosing E c to be 0.2 magnitude units below the maximum plausible magnitude.
Color maps for the maps are adapted from Ahlenius (2005)