Spatial and temporal influence of rainfall on crustal pore pressure based on seismic velocity monitoring

Crustal pore pressure, which controls the activities of earthquakes and volcanoes, varies in response to rainfall. The status of pore pressure can be inferred from observed changes in seismic velocity. In this study, we investigate the response of crustal pore pressure to rainfall in southwestern Japan based on time series of seismic velocity derived from ambient noise seismic interferometry. To consider the heterogeneity of the area, rainfall and seismic velocity obtained at each location were directly compared. We used a band-pass filter to distinguish the rainfall variability from sea level and atmospheric pressure, and then calculated the cross-correlation between rainfall and variations in S-wave velocity (Vs). A mostly negative correlation between rainfall and Vs changes indicates groundwater recharge by rainfall, which increases pore pressure. The correlations differ between locations, where most of the observation stations with clear negative cross-correlations were located in areas of granite. On the other hand, we could not observe clear correlations in steep mountain areas, possibly because water flows through river without percolation. This finding suggests that geographical features contribute to the imprint of rainfall on deep formation pore pressure. We further modelled pore pressure change due to rainfall based on diffusion mechanism. A strong negative correlation between pore pressure estimated from rainfall and Vs indicates that the Vs variations are triggered by pore pressure diffusion in the deep formation. Our modelling results show a spatial variation of diffusion parameter which controls the pore pressure in deep formation. By linking the variations in seismic velocity and crustal pore pressure spatially, this study shows that seismic monitoring may be useful in evaluating earthquake triggering processes or volcanic activity.


Introduction
Pore pressure plays a key role in the occurrence of earthquakes and the volcanic activities (Albino et al. 2018;Ellsworth 2013;Tsuji et al. 2014). Under conditions of critical stress and high pore pressure, small increases in pore pressure can trigger seismicity. Therefore, monitoring the status of pore pressure is a vital part of evaluating dynamic crustal activities. Because pore pressure affects seismic velocity, the state of pore pressure can be assessed by seismic velocity monitoring (Chaves and Schwartz 2016;Hutapea et al. 2020;Nimiya et al. 2017;Rivet et al. 2015;Tsuji et al. 2008;Wang et al. 2017).
In field observations, changes in seismic velocity can be induced by various environmental perturbations (Wang et al. 2017) because seismic velocity is sensitive to variations in stress and water saturation (Grêt et al. 2006). Such perturbations include ocean tides and solid earth tides (Sens-Schönfelder and Eulenfeld 2019), and seismic velocity in coastal locations is sensitive to tidal ocean loading (Yamamura et al. 2003). The influence of the ocean is considered in studies of ambient seismic noise (Hillers et al. 2012). Atmospheric pressure influences seismic velocity over large regions (Niu et al. 2008;Silver et al. 2007), and atmospheric temperature likewise generates seasonal variations in seismic velocity through changes in crustal strain (Ben-Zion and Leary 1986;Berger 1975;Prawirodirdjo et al. 2006), especially in arid regions (Hillers et al. 2015;Richter et al. 2014).
Rainfall and snow are well-known hydrological perturbations by which pore pressure induces seismic velocity changes. For example, the interaction of hydrothermal systems and surface loading from precipitation can lead to seismic velocity reductions (Taira and Brenguier 2016). Snow decreases seismic velocity through increased pore pressure resulting from ice accumulation (Mordret et al. 2016), whereas frost increases seismic velocity at shallow depths by increasing the shear modulus of near-surface materials (Gassenmeier et al. 2015;Tsuji et al., 2012). Rainfall decreases seismic velocity through changes in effective stress (Nakata and Snieder 2012;Miao et al. 2018) and groundwater level (Gassenmeier et al. 2015;Meier et al. 2010;Sens-Schönfelder and Wegler 2006;Tsai 2011). Rainfall triggers seismicity through pore pressure changes caused by crustal loading and unloading (Bettinelli et al. 2008) and pore pressure diffusion (Hainzl et al. 2006;Kraft et al. 2006). Because percolation of water through porous rock may be a contributor to pore pressure changes, we investigated the spatial and temporal relationships between seismic velocity changes and rainfall in a well-instrumented region of Japan.
Crustal deformation in Japan can be evaluated with abundant data from seismic and geodetic observation stations (Aoi et al. 2020). The crust is affected by perturbations from volcanic and seismic activities (Ueda et al. 2013) and surface loads (Heki 2004), including non-tidal ocean loading (Sato et al. 2001). Recent studies have shown that observed seismic velocity changes reflect volcanic activity (Takano et al. 2017;Yukutake et al. 2016) and earthquake activity (Nimiya et al. 2017). Furthermore, seasonal spatial patterns of seismic velocity change throughout Japan can be explained by seasonal variations in rainfall, snow, and sea level (Wang et al. 2017).
This study uses records of seismic velocity changes estimated from ambient noise monitoring in the Chugoku and Shikoku regions of southwest Japan (Fig. 1). This area receives high rainfall from the summer monsoon (Aizen et al. 2001) and is relatively unaffected by volcanic activity and snowfall. To evaluate the influence of rainfall on seismic velocity changes, we performed two-step analyses. In the first step, we sought to identify locations where seismic velocity could be affected by rainfall by directly comparing the seismic velocity to the rainfall via cross-correlations (e.g., Bièvre et al. 2018). The time delay resulting from the cross-correlations helps to constrain near-surface conditions that could be related to lithology-related permeability. In the second step, we modelled pore pressure change due to pore pressure diffusion to estimate a hydrological parameter (i.e. diffusion rate) for the locations where precipitation influence was clearly estimated in the first step. By comparing the seismic velocity change and modelled pore pressure change, we sought to estimate spatial variation of the diffusion parameter in deeper lithology which contributes to predicting precipitation-related pore pressure changes from seismic velocity in Chugoku and Shikoku regions.
Well-quantified monitoring results could be useful information for the evaluation of earthquake triggering mechanisms. In CO 2 geological storage projects and geothermal developments, furthermore, earthquakes induced by fluid injection are a notable public concern. Accurate knowledge of natural pore pressure variations can help in distinguishing whether an earthquake is a natural event triggered by environmental variations or an induced event triggered by fluid invasion.

Data preparation
The Chugoku and Shikoku regions are located in southwest Japan (Fig. 1). The Chugoku region is characterized by mountainous topography with gentle sloping, while the slopes of the mountains in Shikoku Island are mostly steep (Fig. 1b). Figure 1c shows the rock types of our study area from the geological map (Geological Survey of Japan AIST 2015). The Chugoku region is abundant in Cretaceous volcanic and granitic rocks, along with granitic rocks from the Paleocene to the Early Eocene, and Late Pleistocene to Holocene sediments at the northern Chugoku. As for the Shikoku region, Late Cretaceous granite can be found in the northern Shikoku, while Sanbagawa metamorphic rocks are widely distributed across the centre of the Shikoku. Sandstone of Cretaceous-Oligocene accretionary complexes is mainly located in the southern Shikoku Island (steep mountain).
We collected data on seismic velocity changes, precipitation, atmospheric pressure, and sea-level change for the period 2015-2017 in the Chugoku and Shikoku regions. The meteorological data were obtained from the Japan Meteorological Agency (JMA) and we used seismic data from 98 seismometers operated by the National Research Institute for Earth Science and Disaster Resilience (NIED). For each Hi-net station, a three-component sensor of the particle velocity with a natural frequency of 1 Hz is installed in the bottom of the borehole (Obara et al. 2005).
We estimated seismic velocity changes on the basis of ambient-noise coda wave interferometry using the vertical component of ambient noise (Hutapea et al. 2020;Nimiya et al. 2017). To obtain virtual seismograms propagating between pairs of stations, two traces of f A (t) and f B (t) recorded at seismometers A and B were transformed into frequency domain by the Fourier transform: where F A and F B are the seismic waveforms in the frequency domain (ω) recorded at seismometers A and B. The power-normalized cross-correlation (cross-coherence) was applied in the frequency domain between seismometers A and B (e.g., Nakata et al. 2011Nakata et al. , 2015 by where the asterisk (*) denotes a complex conjugate. Changes in seismic velocity between pairs of seismometers were estimated by the stretching interpolation method (Hadziioannou et al. 2009;Hutapea et al. 2020;Minato et al. 2012;Nimiya et al. 2017). This method elongates the time axis and looks for the trace most similar to the reference trace by means of the correlation coefficient CC(ε) between the reference trace and the current trace: where f ref is the reference trace, f cur is the current trace, and t is time. The stretching parameter ε is related to the relative time shift ( �t/t ) and velocity change ( �v/v ) from.
The time window of 100 s for coda waves was used to obtain velocity changes by the stretching interpolation. The seismic velocity change was estimated independently for each individual year by defining the 1-year stack of the coda of cross-correlation data as the reference trace f ref and the 10-day stack of the coda of cross-correlation data as the current trace f cur . To stabilize the monitoring results over the 3-year term, we used the sliding reference method (SRM) to define the reference trace (see Hutapea et al. 2020). In the SRM method, we changed the reference trace for each year. For example, to estimate daily seismic velocity changes in 2015, we defined the coda of cross-correlation stacking over the whole year of 2015 as the reference trace. The daily velocity change was considered to represent the velocity change in the middle of the 10-day window of the current trace. The frequency range of the seismic data was restricted to 0.1 to 0.9 Hz, which reflects the sensitivity of surface waves to S-wave velocity between depths of 1 and 8 km (e.g., Nimiya et al. 2017). To obtain seismic velocity changes for each station (Fig. 1d), we applied spatial averaging within a radius of 40 km.
To obtain precipitation data for each seismic station, we averaged the data from all precipitation gauges within a distance of less than 40 km from the seismic stations (Fig. 1e). Atmospheric pressure and sea-level changes were obtained from the tidal gauge closest to the seismic station ( Fig. 1f ) and the daily sea-level change was estimated by averaging data for the most recent 24-h period.

Methods
Several studies have linked changes in seismic velocity to groundwater recharge by rain precipitation (e.g., Gassenmeier et al. 2015;Sens-Schönfelder and Wegler 2006). When surface water from precipitation replenishes groundwater, we expect decrease in seismic velocity reflecting pore pressure increase due to (a) immediate loading in undrained condition (impermeable), and (b) pore pressure diffusion (Talwani 1997). To confirm the effect of precipitation on changes in seismic velocity, we applied two-step analyses.
In the first step, the time delay between precipitation and seismic velocity change is estimated by cross-correlating the two time series. We identify the locations where velocity changes occurred after precipitations, indicating the influence of pore pressure change due to groundwater recharge. Because the quasi-annual period of seismic velocity change could be influenced by other environmental factors (e.g., atmospheric pressure and sea level), we apply a band-pass filter in order to clearly distinguish rainfall from sea level and atmospheric pressure. Therefore, in the first step, we focus on shorter period fluctuations associated with rain precipitation.
In the second step, we focus on locations where precipitation influence on seismic velocity is clearly identified from the first step. We model pore pressure change based on a diffusion mechanism by groundwater load and compare that with the longer period of seismic velocity change to estimate diffusion rate in deep lithology. Although the observed response is mostly a coupled mechanism (i.e. undrained response and pore pressure diffusion), the effect of pore pressure diffusion may be dominant in the later time, as the pore pressure increase due to diffusion occurs once the immediate loading has dissipated (Talwani 1997). The longer period velocity variation may include the influence of sea level and atmospheric pressure effects, but the seismic velocity variation we used here does not include strong annual features associated with sea level and atmospheric pressure. We summarized the flow of these two-step analyses in the flowchart (Fig. 2).

Step 1: investigation of the rainfall infiltration
To determine an optimal frequency band to clearly distinguish precipitation influences from other environmental factors, we first investigated the power spectra of seismic velocity changes, precipitation, sea-level changes, and atmospheric pressure changes (Additional file 1: Figure S1a). Whereas the power spectrum of seismic velocity changes decreased toward a frequency of 0.1 cycle/day, the spectra of precipitation, sea-level, and atmospheric pressure change showed similar peaks at 0.0018-0.0036 cycle/day, a frequency band close to the annual cycle (Additional file 1: Figure S1b). The similarity of these three peaks meant that the long-term estimated seismic velocity changes could be affected not only by precipitation, but also by sea-level and atmospheric pressure changes.
We excluded frequencies below 0.0036 cycle/day to remove the annual seasonal influence of sea-level and atmospheric pressure changes, and we excluded frequencies above 0.05 cycle/day to eliminate the neap and spring tides of sea-level change and the decreasing spectrum of seismic velocity change. We then searched for the frequency band where precipitation could best be distinguished from sea-level change and atmospheric pressure change, as indicated by weak correlations between precipitation and the other two variables. We applied a band-pass filter for periods between 20 and 137 days (0.05 to 0.0073 cycle/day; Additional file 1: Figures S1c, d) and sought minima in the correlation coefficients between precipitation and sea-level change and between precipitation and atmospheric pressure change, based on the data for all stations.   Figure 3 shows an example of the unfiltered and filtered data for one seismic station (N.YSHH; red dot in Fig. 1f ). Seismic velocity changes in Fig. 3a represent the averaged velocity change within a 40-km radius. The unfiltered data are coloured by the stretching correlation coefficient, averaged for station pairs used to estimate the velocity change. In general, the mean correlation coefficient for station pairs used in this study is above 0.5, and the value is even higher in the period when precipitation is relatively high (e.g., August-November). This indicates the daily seismic velocity change in the study area is stable. Figure 3b-d shows the time series of precipitation, sea level, and atmospheric pressure, respectively. Because the Pearson correlation value between rainfall and sea level is very small (Fig. 3e), we can use band-pass filtering to separate the imprint of precipitation and sea level, as well as precipitation and atmospheric pressure. The correlation coefficients between band-pass filtered precipitation and sea-level and between precipitation and atmospheric pressure change, respectively, are shown in Fig. 4 for all stations. The small Fig. 3 a Example of unfiltered and filtered seismic velocity changes. The colour of the unfiltered velocity change represents the stretching correlation coefficient. b Precipitation, c sea level, and d atmospheric pressure during the study period at the station shown in Fig. 1f. e, f Comparisons of precipitation with changes in sea level and atmospheric pressure, respectively, at the station shown in Fig. 1f. Signals are normalized Fig. 4 Pearson correlation coefficients between band-pass filtered a precipitation and sea-level change and b precipitation and atmospheric pressure change at stations in the study area correlation coefficients indicate that rainfall is distinguishable from sea-level and atmospheric pressure changes.
To further analyse the dependence of seismic velocity changes on rainfall, we applied various time shifts to the rainfall record and evaluated the resulting cross-correlations with seismic velocity changes, as depicted in Fig. 5. Under the assumption that seismic velocity changes are triggered by precipitation after a time lag, we restricted ourselves to positive time lags (i.e. velocity variation after precipitation) and determined the time shift that produced the largest Pearson correlation coefficient.
Although we focus on the shorter cycle (shorter than annual period) in order to clarify the relationship between precipitation and velocity change, it is difficult to distinguish the effects of (a) undrained due to loading and (b) diffusion. Thus, in the next investigation, we evaluate the pore pressure diffusion via modelling.
Step 2: investigation of pore pressure diffusion To calculate the pore pressure change, we use the poroelastic model developed by Talwani et al. (2007). The pore pressure due to diffusion can be described as: where δp k is the water level change, r is the depth from the surface, c indicates the hydraulic diffusion, and δt k indicates the time increment from the starting time k to n , and erfc denotes error complementary function. Although Talwani et al. (2007) also proposed equation for the pore pressure changes due to undrain loading, the effect is smaller than the diffusion in longer period. Here, we consider the contribution of precipitation from 365 days in the past, thus current pore pressure change is calculated by using the summation of pore pressure change from the previous 365 days. We defined water level change from 2015 to 2017 as the deviation from the average precipitation over 2014.
To evaluate the longer period variation, we applied a moving average with 130 days windows for seismic velocity change without band-pass filtering in the first step. By comparing seismic velocity changes with the pore pressure changes computed based on Eq. (7), we estimate the  Fig. 6 Comparison and cross-correlation analysis between seismic velocity changes and precipitation at the station shown as red dot in Fig. 1f: a unfiltered signals, b band-pass filtered signals (normalized), and c band-pass filtered signals with precipitation shifted 8 days later optimum hydraulic diffusion at each station. However, in the calculation of the pore pressure changes, it is difficult to constrain the dependence of the hydraulic diffusion with depth because the relative values of the both parameters in r/ √ c is sensitive to the calculation of the pore pressure by Eq. (7). Therefore, we estimate optimum values of c assuming values of r (i.e. depth), by computing correlation coefficients between observed velocity changes and modelled pore pressure changes. Since we expect decrease in seismic velocity due to increase in pore pressure, we determine the optimum value of c with the largest negative correlation. After optimum values of c are estimated at each station and depth, we construct a map of c.

Results
An example of correlation of rainfall and velocity change at a station near the middle of the study area (red dots in Fig. 1f ) is shown in Fig. 6. Figure 6 includes a correlation without band-pass filtering (Fig. 6a), then after band-pass filtering (Fig. 6b), and then after shifting the precipitation record by 8 days to fit the respective peaks (Fig. 6c). This 8-day time shift raised the correlation coefficient for the 3-year data by more than half, although it is still relatively small at − 0.33. However, when we restricted the comparison to the rainy season (in this case, June to July of 2016 and 2017), the correlation coefficient is much greater (− 0.7).
For most stations there is a negative correlation between rainfall events and seismic velocity changes (Fig. 7a, b); however, a few stations had positive correlations (Fig. 7c). The highest absolute value of these correlations, even after applying the optimum time lag, was approximately 0.3 (Fig. 8). This value is not high because several factors may weaken the correlation between seismic velocity changes and rainfall events. For example, random noise in both of the time series and the time windows decreases the coefficient. In the latter case, a short time window for stacking cross-correlations (10 days in this study) was necessary to analyse the short-term seismic velocity changes induced by precipitation. A longer time window would improve the stability of the velocity change estimate, but would reduce its temporal resolution (Hutapea et al. 2020). Another possibility is that an external factor other than precipitation also influences seismic velocity, such as atmospheric pressure, which can exert effects even at seismogenic depths (Niu et al. 2008). Figure 8a shows the correlation coefficients between rainfall and seismic velocity for all stations, after applying the optimum time lag for each station. Among these stations, 26 stations had absolute correlation coefficients smaller than 0.1, 35 stations had absolute correlation coefficients of 0.1 to 0.2, and 37 stations had absolute correlation coefficients greater than 0.2. We selected the third group with high absolute correlation for further analysis because these stations are clustered and information regarding time lag is reliable only if there is a sufficiently strong correlation. Because most of these stations had a negative correlation between precipitation and seismic velocity, we focused on the stations in this group with negative correlations. These selected stations are shown in Fig. 8b, and their respective time lags are shown in Fig. 8c.
To validate the groundwater recharge due to precipitation, we compared the records of rainfall and groundwater level variations. Figure 9a, b shows the unfiltered records for an example GWL station (see Fig. 1f ) and the calculated cross-correlation between band-pass filtered precipitation and groundwater level. As shown in Fig. 9c, it takes 5 days for rainfall to recharge groundwater, whereas Fig. 9d shows that rainfall is most strongly correlated with a decrease in seismic velocity 9 days later. The small difference in the time lags between Fig. 9c (5 days) and Fig. 9d (9 days) supports our interpretation that the increased groundwater load due to recharge by  Fig. 1c rainfall causes a subsequent decrease in seismic velocity. We show the influence of the near-surface lithology associated with the rainfall infiltration in Fig. 10.
Using the stations where seismic velocity changes are likely influenced by precipitation in Fig. 8, we estimated the optimum value of c by comparing the pore pressure change with seismic velocity change. In Fig. 11, we show an example of comparison between seismic velocity change and pore pressure change with the diffusion rate of 0.14 m 2 /s for depth of 1.5 km that gives the largest negative correlation at the station of N.YSHH (red dot in Fig. 1f ). Figure 12a shows Comparison and cross-correlations between precipitation and ground water level (GWL), and those between precipitation and seismic velocity. The GWL station of this example is shown in Fig. 1f and the seismic station closest to the GWL station is used for comparison. a Relationship between unfiltered precipitation and GWL. b Relationship between band-pass filtered precipitation and GWL. c Relationship between band-pass filtered precipitation and GWL with precipitation shifted earlier by 5 days. d The relationship between band-pass filtered precipitation and seismic velocity change (normalized) with precipitation shifted earlier by 9 days a correlation of the pore pressure and seismic velocity changes for the selected stations in Chugoku and Shikoku regions. Several stations in the Chugoku and Shikoku area show relatively weak negative correlation between pore pressure and seismic velocity change (< 0.2). This can be due to several possibilities; the diffusion is not dominant at these stations, or other perturbations influence the longer-term variations in seismic velocity. Figure 12b shows spatial variation of the estimated diffusion rates for 1-8 km depth, considering the sensitivity depth of surface wave to S-wave in the analysed frequency range (Additional file 1: Figure S3a). The results demonstrate that the hydraulic diffusion controlling the pore pressure spatially varies across the Chugoku and Shikoku regions. The diffusion rates in western Chugoku are generally higher than ones in the eastern area, while station with the highest hydraulic diffusion is found in the eastern Shikoku. Spatial variation of diffusion rate could reflect fracture density. A higher diffusion rate can be interpreted as a well-developed fracture network that connects to a deeper formation.

Influence of near-surface lithology
The delay time of seismic velocity change to rainfall is presumably related to the near-surface conditions, which influence water percolation into geological formation. The time lag between rainfall events and seismic velocity changes in Fig. 8c may represent the time needed for percolating rainfall to reach the water table of an unconfined aquifer. Percolation through the unsaturated zone is likely determined by the permeability of the near-surface layers and the surface geologic and geographic conditions at the seismic station. For example, in mountain regions with high permeability, water derived from the surrounding mountains percolates into intermountain basins. The comparison of our result with the geological and topological map of Japan (Fig. 1b, c) shows that stations with negative correlations are mostly located in granite areas with gentle sloping topography ( Fig. 1b; marked by the colour pink in the legend of Fig. 1c). On the other hand, we cannot identify clear negative correlation in sedimentary rocks with steep slope area in the southern Shikoku ( Fig. 1b; green in Fig. 1c), possibly because water flows away without percolating into deep formation.
Because the unsaturated zone in humid climates is generally less than 10-m thick (Phillips and Castro 2003), we assumed the unsaturated zone in our humid study area to be shallower than 10 m. Borehole logs from the sites where our seismometers are deployed classify the shallow formation as high-permeability materials and weathered igneous rocks (Obara et al. 2005). Under the assumption that S-wave velocity may be related to permeability, we examined plots of time lag versus S-wave velocity (Fig. 10) to evaluate the relationship between lithology and time lag. Although the relationships are unclear, we identified some features for each formation.
Among the 29 stations obtained from the step 1, the lithology of 19 stations can be classified into highpermeable materials (Fig. 10a) and weathered igneous rocks (Fig. 10b). A total of 8 stations with high-permeability material such as sandy soil, silt, and gravel shows a positive trend in which the time lag increases with increasing S-wave velocity (Fig. 10a). Because the seismic velocity varies inversely with porosity, this relationship confirms that percolation could be faster in more porous materials (i.e. low Vs) and slower where porosity is lower.
In weathered igneous rocks, which consist mostly of granite, the time lag and S-wave velocity show a modest negative trend (Fig. 10b). This trend may be connected to the spatial concentration of fractures in these rocks. Although fractures in crystalline rocks generally decrease with increasing depth, fractures are the primary determinant of permeability at depths shallower than 10 m in plutonic and crystalline metamorphic rocks (Freeze and Cherry 1979). Furthermore, the decreasing time lag with increasing S-wave velocity implies that the water table is shallower in the less-fractured igneous rocks, whereas brittle, morefractured igneous rocks allow rainfall to percolate to greater depths, resulting in longer lag times for water to reach the saturated zone.
The estimated time lag can be also influenced by the delayed response of pore pressure change associated with the diffusion mechanism. Indeed, the time lag between seismic velocity change and rainfall is longer than that between groundwater level and rainfall (Fig. 9). This might reflect the influence of the delayed response of pore pressure change (i.e. seismic velocity change), in addition to the time delay due to percolation of rainfall to the water table. Fig. 11 Comparison of moving averaged seismic velocity changes and pore pressure estimated from precipitation at the seismic station (red dots in Fig. 1f ). a Moving averaged seismic velocity change. b Precipitation and the calculated pore pressure change. c Correlation between averaged seismic velocity change and pore pressure. The signals are normalized Fig. 12 a The correlation map between seismic velocity change and pore pressure. b The map of diffusion parameter for each depth (1, 1.5, 2, 4, 6, and 8 km). The stations with negative correlations smaller than 0.2 are not included on the map. The colour bar in each panel represents the different range of hydraulic diffusion rate Seismic velocity changes due to pore pressure diffusion The example shown in Fig. 11 demonstrates that pore pressure increases from July to November. This pattern agrees with the seismic velocity decrease from July to November. A similar pore pressure and seismic velocity variation also occurs at other locations across Chugoku and Shikoku (Additional file 1: Figure S3). It is known that rainfall in July-September can trigger seasonal seismicity in the Chugoku area (Ueda and Kato 2019). The similar timeline suggests that the longer period seismic velocity change might have been influenced by pore pressure change induced by rainfall, although the long period variation is also influenced by sea level and atmospheric pressure variations.
The surface wave depth sensitivity to S-wave could be associated with the frequency of the coda wave used to estimate seismic velocity change (Nimiya et al. 2017). Although the frequency range of our seismic velocity change is sensitive to 1-8 km, the largest sensitivity derived from velocity model in our study area seems to be within 1.5-2 km depth (see Additional file 1: Figure S2). Within these depths, the range of hydraulic diffusivity for Chugoku and Shikoku regions varies from 0.02-1 m 2 /s (Fig. 12b). Suppose we take an example of 1.5 km depth, then the hydraulic diffusion rate at the western Chugoku area would be 0.09-0.14 m 2 /s, 0.09-0.25 m 2 /s for northern Chugoku, 0.02-0.1 m 2 /s for eastern Chugoku, and 0.02-0.5 m 2 /s for eastern Shikoku.

Mechanisms of pore pressure variation
We summarize our results and interpretation in Fig. 13. The near-surface condition (e.g., lithology and fracture) controls the percolation of rainfall, resulting in a time lag between rainfall and pore pressure increase (Fig. 13a). After the rain precipitation, we expect pore pressure increase mainly due to (a) immediate loading in undrained condition and (b) pore pressure diffusion.
As groundwater level increases due to rainfall infiltration, immediate loading causes pore pressure increase from Pp1 to Pp2 in Fig. 13a, and generate thin cracks (white arrow in Fig. 13b). This condition persists until pore pressure ceases to the surrounding fractures in deep formation (Pp2 to Pp3 in Fig. 13a). This pore pressure variation could be mainly observed by shorter period seismic velocity reduction (Fig. 6). Then, the load from groundwater level increase triggers pore pressure diffusion through the pre-existing fracture network. As the pore pressure front arrives (white wavy arrow in Fig. 13c), there is an increase of pore pressure from Pp4 to Pp5 in Fig. 13a. This pore pressure increase can be monitored by longer period seismic velocity (Fig. 11).
We conclude that local lithology, both above the groundwater table and in the deep formation, contributes to the pore pressure changes associated with rainfall. The interpretations we describe here are simple ones. In real hydrogeological systems, however, there are many other complex mechanisms that affect the time lag (e.g., flow path influenced by geographical features), as well as the fracture permeability in the deeper lithology.

Conclusion
The status of pore pressure changes associated with rainfall can be evaluated by monitoring the seismic velocity. By calculating the cross-correlation between rainfall and seismic velocity changes, we can identify the locations where seismic velocity change is influenced by precipitation. Furthermore, by modelling pore pressure change based on pore pressure diffusion due to rainfall, we can constrain hydraulic diffusion from long-period seismic velocity changes. Our primary conclusions are: 1. The influence of rainfall on seismic velocity change varies depending on the lithology. Clear negative correlations between rainfall and seismic velocity can be observed in the granite areas and terrains with gentle topography. On the contrary, there are no clear correlations observed in the steep mountain areas. 2. The time lag between precipitation and seismic velocity change constrains near-surface conditions that could be related to lithology-related permeability. Similar time lag between precipitation and (See figure on next page.) Fig. 13 Summary of the mechanism of crustal pore pressure change (Pp) associated with rainfall (modified from Talwani et al. 1997). a The time duration for the immediate high pore pressure (due to undrained effect) and pore pressure diffusion to occur with the increasing water level. The schematic figures of b immediate loading and c pore pressure diffusion in a later time. The white arrow in b represents the immediate increase of pore pressure due to undrained condition. The white wavy arrows in c represent the pore pressure diffusion ground water level demonstrates that the increased groundwater load causes a subsequent decrease in seismic velocity. 3. The pore pressure diffusion caused by rainfall infiltration can be modelled and controls longer term pore pressure change. The spatial variation of diffusion parameter estimated by the modelling depends on fracture connectivity and is spatially varied.
Additional file 1: Figure S1. (a) Map of seismic stations. (b) The power spectra of seismic velocity changes, precipitation events, sea-level changes, and atmospheric pressure changes for the 0-0.5 cycle/day frequency band for the seismic station shown in red of panel (a). Precipitation is averaged around the seismic station, and sea level and atmospheric pressure are taken from the closest sea-level and pressure gauges. (c) The shape of weighting functions of the band-pass filter for selected frequencies defined between four points (f1 = 0.0073, f2 = 0.0082, f3 = 0.03, and f4 = 0.05 cycle/day). (d) The band-pass filter applied to the power spectra. The unshaded parts of the power spectra were used in the analysis. The peaks with the lowest frequency represent annual or quasi-annual cycles. Figure S2. Depth sensitivity of surface wave (Rayleigh wave) to S-wave velocity. The surface wave sensitivity to S-wave was calculated by DISPER80 (Saito 1998) for 1D velocity layer model at the Chugoku region (Nishida et al. 2008). The amplitudes were normalized by the maximum sensitivity of frequency 0.9 Hz. Figure S3. Comparison of moving averaged seismic velocity changes and the calculated pore pressure at the seismic stations of A and B. The top panel indicates the moving averaged seismic velocity change. The middle panel shows the comparison of precipitation and pore pressure. The bottom panel represents the correlation coefficient between averaged seismic velocity change and pore pressure. The signals are normalized.