Seasonal Modulation of Crustal Seismicity in Northeastern Japan Driven by Snow Load

Numerous studies have reported that surface hydrological loading can seasonally modulate seismicity rates at crustal depths. For example, substantial winter snow accumulation occurs across the Japanese Islands, and these snowy regions appear to have seasonally modulated the occurrence of previous large inland earthquakes. Therefore, it is important to investigate the impact of seasonal stress changes on crustal seismicity to deepen our understanding of earthquake generation. Here we constrain seasonal changes in the surface load across northeastern Japan using Global Navigation Satellite System surface displacements and evaluate the potential relationship between temporal trends in inland seismicity and estimated seasonal stress changes. The spatial distribution of the seasonal surface load is consistent with snow depth along the Sea of Japan. The inland seismicity beneath northeastern Japan is modestly modulated by the seasonal stress changes that are induced by the annual snow load. However, this seasonal response is weaker than that in other regions. This weak modulation may be due to the small surface‐load‐induced stress perturbation relative to the long‐term‐averaged stressing rate and/or the limited presence of crustal fluids to trigger seismicity in Japan.


Introduction
Seismicity rates may correlate with natural phenomena at the surface/near surface that produce temporal changes in the local stress field and fault strength, such as hydrological loading (Amos et al., 2014;Bettinelli et al., 2008;Craig et al., 2017;Heki, 2003;Hsu et al., 2021;Johnson et al., 2017bJohnson et al., , 2017aJohnson et al., , 2020;;She et al., 2022;Xue et al., 2020Xue et al., , 2021;;Yao et al., 2022) , precipitation (Hainzl et al., 2006(Hainzl et al., , 2013)), atmospheric pressure changes (Gao et al., 2000), and surface temperature variations (Ben-Zion & Allam, 2013).For example, the seismicity rate in California varies seasonally in response to stress changes related to annual variations in hydrological loading (Amos et al., 2014;Johnson et al., 2017b).Johnson et al. (2017b) calculated the stress changes that hydrological loading imposes on fault planes based on Global Positioning System (GPS) vertical displacements and showed that the seismicity rate in California is more likely to increase when the shear stress increases.Seasonal variations in seismicity rates have also been observed in Japan, and the seasonal snow load and precipitation rate have been suggested as plausible causes (Heki, 2003;Ohtake & Nakahara, 1999;Okada, 1982;Ueda & Kato, 2019;Xue et al., 2021).Heki (2003) reported that previous large inland earthquakes (𝑀 ≥ 7) in the snowy regions of the Japanese Islands tended to occur more frequently in spring and summer than in fall and winter owing to the spring thaw potentially inducing the unclamping of fault planes.However, few studies have directly compared the seasonal variations in crustal seismicity in Japan with surface-load-induced stress changes.
Heavy snow accumulates on the western flanks of the mountain ranges in northeastern Japan during the winter months, and this snow loading generates annual displacements that are observable at local GPS sites (Heki, 2001).Here, we estimate the monthly seasonal surface load using long-term Global Navigation Satellite System (GNSS) surface displacements and comprehensively evaluate the temporal trends in inland seismicity with the estimated seasonal stress perturbations that are resolved on the nodal planes of typical fault mechanisms.

GNSS and snow data
We use the final solutions of the GEONET (GNSS Earth Observation Network) daily coordinates, which are also known as the F5 solutions (Takamatsu et al., 2023), that were determined for the 2004-2010 period (Figure 1).We only use the vertical displacements, which are more sensitive to the surface load than the horizontal displacements, to estimate the seasonal surface load.

Maximum snow depth estimates from the AMeDAS (Automatic Meteorological Data
Acquisition System) sites determined for the 2004-2010 period are used for comparison with the estimated seasonal surface load.

Vertical GNSS time series and surface-load modeling
The raw GEONET time series contains common-mode noise at all stations, and this noise may obscure the appropriate seasonal displacement variations.We remove the common-mode errors (CMEs) using a procedure that employs principal component analysis (PCA)-based spatiotemporal filtering (Dong et al., 2006).First, we remove the sudden displacements caused by   ≥ 6.0 (  ; Japan Meteorological Agency (JMA) magnitude) earthquakes that occur within a distance  and are expected to possess a displacement  of 1 cm or more (Okada, 1995): (1) Furthermore, we remove any sudden displacements that occurred during antenna replacements (dates are published by GSI).We adopt PCA to a matrix , whereby each column contains the detrended and demeaned vertical displacements from a single station (  ), and the rows contain the displacements for all stations at a given epoch (  ).The matrix  is decomposed by: where   () is the th principal component, which represents the temporal variations, and   () contains the eigenvectors, which represent the corresponding spatial responses to the principal components.We define the first principal component, which is the most spatially uniform mode, weighted by the medium of  1 () as the CME.
We use the GNSS stations in southern Japan (cyan triangles in Figure S1), where it rarely snows, to estimate CME.The estimated CME time series is shown in Figure 1b, and example raw and CME-corrected displacement time series are shown in Figures 1c and d, respectively.
The CME-corrected data are used to extract seasonal variations via the Greedy Automatic Signal Decomposition algorithm (Figure 2) (Bedford & Bevis, 2018).We assume that the verticaldisplacement time series can be expressed in the following form:  vertical displacements, which were estimated using stations distributed across southern Japan (cyan triangles in Figure S1).Time series of (c) raw and (d) CME-corrected vertical displacement at station 970810 (red triangle in Figure 1a).
We select the seasonal vertical displacements (e.g., Figure 2d) extracted from 303 stations (yellow triangles in Figure 1a) to estimate the seasonal surface load (monthly averages) during 2004-2010.We exclude the vertical displacements from ten stations in this analysis because the displacements at five stations are closely related to groundwater-level changes due to local irrigation (blue triangles in Figure 1a; Munekane et al., 2004;Tobita et al., 2004) and the displacements at another five stations contain high seasonal amplitudes (green triangles in Figure 1a).The extracted seasonal components of the vertical displacements are inverted to obtain the seasonal surface loads as a function of location (Figure S2) by simultaneously minimizing the misfit of the model to the data, a Laplacian term that constrains the spatial variations in the surface load between neighboring grid points, and a dumping term that limits the surface-load changes.The surface loads are estimated at 0.25° intervals (latitude and longitude).The Green's function of the vertical displacements   (m) at Earth's surface represents a circular disc load (1 kg/m 2 ) with an angular radius  (°), which takes the following form (Wahr et al., 2013): where: ℎ  contains the load Love numbers, which are calculated using the approach in Wang et al.
(2012) and the Preliminary Reference Earth Model (PREM; Dziewonski & Anderson, 1981),  is Newton's gravitational constant,  is Earth's radius,  is the gravitational acceleration at Earth's surface,  is the angular distance away from the center of the disc load, and   is the Legendre function.We set the angular radius  to 0.125°, then minimize the misfit function as follows: where  is the data vector, which contains the extracted seasonal component from the observed GNSS vertical displacements,  is a weighting matrix (  =     ⁄ , where   is the root-meansquare amplitude of the residual term of each station, Equation 3),  is the model vector, which contains the model surface load at each grid point,  is the matrix consisting of Green's functions (Equation 4), and  is the Laplacian operator.We employ a grid search scheme to the estimated  1 and  2 values that minimize Akaike's Bayesian Information Criterion (ABIC; Akaike, 1980), following Fukahata et al. (2004).

Surface-load results
The spatial distribution of the seasonal changes in the March-August vertical displacements is shown in Figure 3a.The spatial distribution of the surface load in Figure 3d is estimated using these displacements (Figure 3a), with ~0.5° (latitude and longitude) spatial smoothing applied.A positive load is estimated in areas of subsidence, and a negative load is estimated in areas of uplift.Furthermore, areas with a positive load (Figure 3d) are spatially coincident with the snowy region (Figure 3f).
The dependence of  1 and  2 in Equation 6on ABIC of the surface-load models is shown in Figure S3.The difference in ABIC from the best (smallest ABIC) model (ΔABIC) indicates that exp(−∆ABIC/2) can be interpreted as the relative probability of a given model fit compared with the best model (Akaike, 1980;Kumazawa & Ogata, 2013).Therefore, models with ΔABIC < 6 can be accepted with a probability of >5%.The spatially roughest and smoothest surfaceload models with ΔABIC < 6 are shown in Figure S4a and b, respectively.The spatial distributions of the two models are similar to that of the best model (Figure 3d), which suggests that the best model provides a robust estimate of the surface-load distribution.0.98 kPa.The estimated surface load fluctuates around zero from June to October (average of 0.01 ± 0.66 kPa; Figures 4f-j), when the region is snow free.The opposite trend is generally observed in the monthly estimated surface loads along the Pacific Ocean, but with lower amplitudes.This trend is inconsistent with the increase in surface load due to snowfall during the winter season.Our results reveal that the positive surface loads (Figure 3d) correspond to the snowy region in northeastern Japan (Figure 3f).We demonstrate that the snow load is a plausible source of the vertical subsidence by comparing the spatial distribution of the estimated surface load (Figure 3d), which was derived from the observed vertical displacements (Figure 3a), with the spatial distribution of the snow load by spatially smoothing the observed averaged maximum snow depth data for March at the AMeDAS sites (Figure 3f).There is a positive correlation between the estimated surface load and observed snow depth when a snow density of 400 kg/m 3 is assumed, as shown in Figure 3e.
We also calculate the expected vertical displacement at each GNSS station (Figure 3c) based on the Green's functions (Equation 4), whereby we use the averaged maximum snow depth (Figure 3f) and an assumed snow density of 400 kg/m 3 to calculate the observed snow load.The spatial distribution of the expected vertical displacement (Figure 3c) roughly corresponds to the amplitude of the observed subsidence at the GNSS stations (Figure 3a) along the Sea of Japan.However, the maximum vertical subsidence (~6 mm) that is expected from the snow depth is less than the maximum observed GNSS vertical displacement (~15 mm).This suggests that the observed GNSS vertical subsidence includes loading sources other than the snow mass, which is obtained from the snow depth.For example, geodetic investigations of groundwater from snow melt along roads in the Niigata area have indicated significant subsidence (Morishita et al., 2020;Sato et al., 2003;Shimada et al., 2021).Or the displacement measurements may be incorrect because of the accumulation of snow on the GNSS pillar and/or radome covering the antenna, thereby impeding the GNSS signals and enhancing signal scattering effects (Jaldehag et al., 1996).Another source of this uncertainty may be the regularization of the current inversion scheme, which smooths the surface-load model and makes it less sensitive to the data variability.Although factors other than snow load may contribute to the observed vertical displacement, we conclude that snow load is a plausible cause of the seasonal variations in vertical displacements along the eastern margin of the Sea of Japan.
The seasonal component of the vertical displacement along the Pacific Ocean from March to August shows signal consistent with uplift (Figure 3a), and a corresponding negative surface load is estimated (Figure 3d).Although this result is inconsistent with snow accumulation in winter, the uplift in March has been confirmed by Fujiwara et al. (2022), who extracted the averaged displacement in March after the 2011 Tohoku-Oki earthquake.The observed uplift along the Pacific Ocean is considered a true signal; however, it is uncertain whether this signal can be interpreted as deformation owing to changes in the surface load.We discuss the contribution of the changes in surface load along the Pacific Ocean margin in Section 3.4 to evaluate the relationship between surface-load-induced stress changes and temporal trends in inland seismicity.
In this study, we extract the averaged seasonal component during the 2004-2010 period (Figure 2d).We observe that the magnitude of the vertical subsidence increases as more snow accumulates across the region (Figure S5).We note that any seasonal changes in the timevarying amplitudes (Cleveland et al., 1990;Köhne et al., 2023) should be examined in a future study.Furthermore, careful consideration must be taken to exclude the effects of both groundwater (Morishita et al., 2020;Sato et al., 2003;Shimada et al., 2021) and snow accumulation on the GPS antenna (Larson, 2013;Larson et al., 2015) from the vertical displacements when estimating the spatiotemporal distribution of the surface load.
3 Evaluation of temporal trends in inland seismicity and resolved stress perturbations

Earthquake catalog data
The earthquake data in our analysis are taken from the JMA hypocenter catalog (Figure 5).We use   ≥ 3.0 shallow crustal earthquakes (≤25 km depth) that occurred during 1980-2010.We confirm that the   ≥ 3.0 earthquakes are almost completely detected based on the obtained goodness of fit to Gutenburg-Richter law of >90 (Figure 5b) (Wiemer & Wyss, 2000).
Aftershocks should be removed via a declustering technique prior to investigating the seasonal modulation of seismicity.Here we apply the Hierarchical Space-Time Epidemic-Type Aftershock Sequence (HIST-ETAS) model (Ogata, 2004;Ogata et al., 2003;Ueda et al., 2021) to the earthquake catalog and the probability that each event is a background event is used as the declustered catalog (Ueda & Kato, 2023;Zhuang et al., 2002).A comparison of the cumulative numbers of total and background earthquakes (Figure 5c) shows that the aftershocks are successfully removed from the original earthquake catalog by the HIST-ETAS model.

Seasonal modulation of inland seismicity
We focus on the spatial distribution of the monthly surface load (Figure 4), which was derived from the vertical displacements (Section 2.2), to evaluate the relationship between the seasonal surface-load-induced stress changes and the temporal trends in inland seismicity.We employ the approach in Johnson et al (2017b), whereby we calculate the strain rate tensors using a modified version of STATIC 1D (Pollitz, 1996) that is adapted for a vertical force at the surface (Pollitz et al., 2013) and an assumed PREM structure (Dziewonski & Anderson, 1981).
The strain rate tensors are calculated at 10 km depth, which is considered a representative seismogenic depth for the study area.We then convert the strain rate tensors to stress rate tensors by assuming isotropic elasticity with a Poisson ratio of 0.25 and a shear modulus of 35 GPa.We calculate the shear, normal, and Coulomb stressing rate (for friction coefficients  = 0.1, 0.4, and 0.7) on the two planes of maximum shear stress (the two nodal planes where the shear stress is at a maximum in the tectonic stress field) (Figure 6; Terakawa & Matsu'ura, 2008, 2010) at 0.05° intervals (latitude and longitude).We compute the seasonal stress changes by integrating the calculated seasonal stressing rates.The spatial distribution of the monthly Coulomb stressing rates and monthly Coulomb stress changes ( = 0.4) are shown in Figures 7 and 8, respectively.
We evaluate the temporal trends in inland seismicity with respect to the seasonal surfaceload-induced stressing rate/stress changes based on the percentage of excess earthquakes that occurred during a range of stress intervals (Cochran et al., 2004;Johnson et al., 2017bJohnson et al., , 2017aJohnson et al., , 2020;;She et al., 2022;Xue et al., 2021).The percentage of excess earthquakes is calculated according to the following equation: where   is the percentage of excess earthquakes per stress interval,   is the actual number of events per stress interval, and   is the expected number of earthquakes under the assumption of a temporally uniform distribution.The number of earthquakes beneath the region Figure 6.Spatial distribution of planes of maximum shear stress in the tectonic stress field (Terakawa andMatsu'ura, 2008, 2010).
where we estimated the surface load is weighted by the probability that each event is a background event to remove the effect of aftershocks (Figure 5a).We assume that each earthquake is driven by either of the two planes of maximum shear stress at the nearest calculated grid point.We obtain   by first randomly selecting one of the two nodal planes (Figure 6) and then adopting the values of stress change or stressing rate in the month when an are binned in the outermost stress and stressing rate intervals.We calculate   (Equation 7) and its standard deviation by making a total of 500   and   samples.

Stress modulation of inland seismicity
The relationships between excess seismicity and stressing rates and stresses induced by 266 and stress changes (Figure 10) is more stable than that between excess seismicity and stressing rates (e.g., weak negative correlation with the shear stressing rate (Figure 9b)).The slopes of the percent excess seismicity to the stressing rate and stress are ~0.4N Ex /(kPa/yr) and ~3 N Ex /kPa, respectively.The relationship between the percent excess seismicity and seasonal surface-loadinduced stressing rate/stress has been explored in California (Johnson et al., 2017a, b), Alaska (Johnson et al., 2020), and the Lake Biwa region (Xue et al., 2021).The slopes of the percent excess seismicity to the stressing rate and stress in the inland Tohoku region are smaller than or comparable to those in other studies (2-5 N Ex /(kPa/yr) (Xue et al., 2021), >3 N Ex /kPa (Johnson et al., 2017b)), which suggests that the seismicity in the inland Tohoku region is less sensitive to the surface-load-induced stress changes compared with that in other regions.

Discussion
We calculate the seasonal stress changes using the surface load (Figure 4), which is estimated from the seasonal component of the vertical displacement.However, it is unclear whether the seasonal variations in vertical displacements along the Pacific side of northeastern Japan reflect surface-load changes because they are inconsistent with the snow-load trend, as mentioned in Sections 2.3 and 2.4.Therefore, we evaluate the relationships between the percent excess seismicity and calculated stressing rate and stress changes using the spatial distribution of the surface load along only the Sea of Japan side of northeastern Japan (Figures S6 and S7,respectively).We do not identify a correlation between the percent excess seismicity and stressing rates/stress changes beyond the standard deviation, which suggests that the surface-load changes on the Pacific side may be important in modulating inland seismicity.The mechanism(s) for these surface-load changes should be considered in future studies.
Here we assume that the seasonal surface load is the primary factor that modulates the seasonal variations in inland seismicity and investigate the relationship between inland seismicity and surface-load-induced stresses.However, there are other possible mechanisms, such as heavy precipitation infiltration, atmospheric pressure changes, and surface temperature variations (e.g., Ben-Zion & Allam, 2013;Hainzl et al., 2013;Johnson et al., 2017b).Therefore, we also evaluate the effect of ocean loading, which is excited by interactions between surface winds, atmospheric pressure, temperature and density gradients, and currents (Dong et al., 2002).
We use the monthly averaged dynamic sea surface height data calculated by the Estimating the Circulation and Climate of the Ocean Project (http://www.ecco-group.org);data are available through the Jet Propulsion Laboratory (https://podaac.jpl.nasa.gov).The dynamic sea surface height around Japan is low in spring and high in fall.We convert the monthly changes in the dynamic sea surface height to ocean mass loads by multiplying the heights by the density of water (1000 kg/m 3 ).The ocean mass load is discretized on a 0.25° grid (latitude and longitude), and the stressing rate and stress change on the nodal planes induced by the seasonal ocean mass load are calculated using the method in Section 3.2.The amplitude of the seasonal variations in the ocean mass load-induced stress changes is several tens of Pa, which is quite small compared with that caused by the terrestrial surface load (several kPa, Figure 8).The relationships between the percent excess seismicity and stressing rates and stresses caused by the summation of terrestrial surface load (Figure 4) and ocean mass load are shown in Figures S8 and S9, respectively.The positive correlations between the seismicity and stressing rates/stress are still visible, even when the effect of the ocean mass load is considered.(Johnson et al., 2017b).Therefore, the hydrological load is likely to be the largest source of seasonal stress, even in Japan; however, additional research is required to verify the importance of the hydrological load on the seasonal stresses in Japan.
We find that the seismicity (M ≥ 3.0) in the inland Tohoku region is modestly modulated by the surface-load-induced seasonal stressing rate/stresses (Figures 7 and 8).Heki (2003) reported that the historical large earthquakes (M ≥ 7.0) in snowy region occurred more frequently in spring and summer than in fall and winter, whereas the bimonthly distribution of moderate earthquakes (6.0 ≤ M ≤ 7.0) in the snowy region did not show a clear seasonal variation.However, the location accuracies of the epicenters of historical earthquakes along the Sea of Japan are low, such that these earthquakes may not have truly occurred beneath the snowy region owing to limited information on the associated seismic damage in historical documents.
Furthermore, Heki (2003) did not quantitatively evaluate the relationship between the temporal trends in inland seismicity and snow-load-induced stress changes.We identify a positive correlation between the excess seismicity and snow-load-induced stress changes, which supports the idea that large earthquakes are also seasonally modulated by snow load (Heki, 2003).
Therefore, the seasonal variability in seismicity in the inland Tohoku region may be a scaleinvariant phenomenon.This phenomenon has also been observed in San-in district, southwest Japan (Ueda & Kato, 2019) and California (Johnson et al., 2017a).
Our results (Figures 9 and 10) suggest that seismicity in the inland Tohoku region is less sensitive to surface-load-induced stress changes than that in other regions (Johnson et al., 2017a(Johnson et al., , 2017b(Johnson et al., , 2020;;Xue et al., 2021).This reduced sensitivity may indicate that the long-termaveraged stressing rate in the inland Tohoku region is higher and the stress change perturbation (ratio of the seasonal stressing rate to the long-term-averaged stressing rate) is smaller than that in other regions, which is supported by the presence of strain concentration zones in the Niigata-Kobe Tectonic Zone and Ou Backbone Range (Miura et al., 2002(Miura et al., , 2004;;Sagiya et al., 2000).
Another controlling factor may be that there are fewer crustal fluids to trigger seismicity in Japan (Brodsky et al., 2003;Harrington & Brodsky, 2006).Harrington and Brodsky (2006) found that the seismicity in Japan is less dynamically triggered by distant earthquakes than that in other regions; one possible mechanism for this reduced dynamic triggering is that the seismic waves from teleseismic earthquakes cannot unclog fluid-saturated crustal fractures, resulting in the relative inability of fluid infiltration into the surrounding pre-existing faults.The internal fluid pressure on the faults cannot increase because of frequent fracturing from nearby earthquakes before the fluids can saturate the cracks.This mechanism may be applicable where there are seasonal variations in seismicity if the fractures can be unclogged by the reduced surface loading, as opposed to the propagating seismic waves from teleseismic earthquakes.

Conclusions
We estimated the spatial distribution of the surface load using GNSS vertical displacement data from the inland Tohoku region and examined the relationship between subsurface stress changes induced by the estimated seasonal surface load and temporal trends in inland seismicity.The seasonal surface load increases in winter, decreases in spring, and is almost zero in summer and fall along the Sea of Japan.This seasonal surface load trend corresponds to the observed snow depth changes at the AMeDAS sites.Although we find a positive correlation between excess seismicity and stress changes, the seasonal modulation of seismicity is weaker than that observed in other regions.This might be due to the smaller perturbation of the stress changes relative to the long-term-averaged stressing rate and a relative inability to unclog the fluids in crustal fractures by reducing the surface load.
This study was supported by JSPS KEKENHI Grant Numbers 20J11654 and 22KJ1770.This research is a part of the PhD thesis of Taku Ueda (Ueda 2022).
Figure1d) without assuming the timing and number of any step or transient changes(Bedford & Bevis, 2018).

Figure 1 .
Figure 1.GNSS station distribution and example vertical-displacement time series.(a) Map of the GEONET GNSS station distribution used in this study.Yellow triangles denote the GNSS stations used to estimate seasonal surface loads.The red triangle is station 970810, which is a representative example used in Figures 1, 2, and S5.Blue triangles denote stations whose displacement signals were closely related to groundwater-level changes due to local irrigation.Green triangles denote stations whose displacement signals suffered from extremely high seasonal amplitudes.(b) Time series of common-mode errors (CMEs) for vertical displacements, which were estimated using stations distributed across southern Japan (cyan triangles in FigureS1).Time series of (c) raw and (d) CME-corrected vertical displacement at station 970810 (red triangle in Figure1a).

Figure 2 .
Figure 2. Decomposition of the vertical-displacement time series.(a) CME-corrected vertical displacements at station 970810 (red triangle in Figure 1; same figure as Figure 1d).(b) Trend and transient component, (c) step component, (d) seasonal component, and (e) residual of the CME-corrected vertical displacement in (a), with each component extracted via the approach in Bedford and Bevis (2018).

Figure 3 .
Figure 3. Spatial distributions of vertical displacements, surface load, and maximum snow depth across the study area.(a) Extracted seasonal component of the vertical displacements in mid-March (relative to mid-August).(b) Calculated vertical displacements using the spatial distribution of surface load in (d).(c) Expected vertical displacements using the observed snow load, which was derived from the averaged maximum snow depth in (f).(d) Spatial distribution of surface load, which was estimated from the displacements in (a).(e) Relationship between the estimated surface load in (d) and averaged maximum snow depth in (f) for every grid point in the model.The dashed line shows the theoretical relationship between snow depth and surface load using a snow density of 400 kg/m 3 .(f) Averaged maximum snow depth data for March based on the AMeDAS site observations (black squares).

Figure 4 .
Figure 4. Spatial distribution of monthly surface load, estimated using the seasonal component extracted from the vertical displacements.The black line marks the boundary between the area along the Sea of Japan and the area along the Pacific Ocean.

Figure 5 .
Figure 5. Earthquake data used in this study.(a) Epicenter distribution (  ≥ 3.0).Colors indicate the probability a given event is a background earthquake by applying the HIST-ETAS model (e.g., Ogata 2004).(b) Cumulative frequency-magnitude distribution.(c) Magnitude-time plot (green circles).Red and blue lines denote the cumulative numbers of total earthquakes and background earthquakes, respectively.

Figure 7 .Figure 8 .
Figure 7. Spatial distributions of the monthly Coulomb stressing rates ( = 0.4) on the planes of maximum shear stress at 10 km depth induced by the seasonal surface load.

Figure 9 .
Figure 9. Relationships between excess seismicity and stressing rates.Percent excess seismicity for the (a) normal stressing rate, (b) shear stressing rate, and Coulomb stressing rates for (c)  = 0.1, (d)  = 0.4, and (e)  = 0.7.The bin interval is 6 kPa/year.Error bars show 1σ standard deviation.The slope of the best-fit line and its 1σ standard deviation are shown above each panel.

Figure 10 .
Figure 10.Relationships between the excess seismicity and stresses.Percent excess seismicity for the (a) normal stress, (b) shear stress, and Coulomb stresses for (c)  = 0.1, (d)  = 0.4, and (e)  = 0.7.The bin interval is 1 kPa.The error bars show the 1σ standard deviation.The slope of the best-fit line and its 1σ standard deviation are shown above each panel.
Johnson et al. (2017b)evaluated the stresses caused by various types of loading sources, such as hydrological loads, atmospheric pressure changes, temperature variations, and Earth pole tides, and concluded that hydrological loads are the largest source of seasonal stresses in California.The stresses induced by the seasonal surface load and ocean mass load in the present study are comparable to those induced by the estimated hydrological load and ocean mass load in California