Abnormal Strain Induced by Heavy Rainfall on Borehole Strainmeters Observed in Taiwan

: We found some obvious abnormal strain induced by heavy rainfall from borehole strainmeters deployed in Western Taiwan. The strain induced by rainfall can be divided into two parts, one is the quick response for extra loads of rainwater on the ground, and another one is the slow response for rainwater inﬁltrating into the strata. The quick and slow rainfall responses of areal strain data are analyzed using the technique of recursive digital ﬁltering. Moreover, the rainfall impact functions of the studied stations are calculated using deconvolution. We found, in most cases, the response strain will reach the maximum in half an hour after heavy rainfall, and then show an exponential decay, it might persist more than 200 h depending on the hydrogeological condition around the station. Whereas the river ﬂowing beside the station will help accelerating the runoff dispersion and reducing rainfall decay time in the hill or mountain region. We also compare the results after calibration in term of isotropic and vertical coupling individually. We found that the response strains are smaller in vertical coupling rather than isotropic coupling. The effects of debris avalanches caused by intensive rainfall in the mountain areas can be viewed as two types of rock deformation: generated only under the inﬂuence of rainfall and generated by the increased load in the river channels due to rainfall-induced landslides or debris ﬂow. When the cumulative rainfall exceeds a certain threshold, the strain response curves show a noticeable anomaly likely due to the effects of the debris ﬂow events in places prone to landslides.


Introduction
Strainmeters feature high precision and wide frequency ranges and are able to record the crustal deformation caused by tectonic activities. Many of them have been installed in seismically active tectonic settings around the world such as America, Japan and Taiwan. Apart from the application of observations and research of seismic precursors, strainmeters are also frequently utilized to observe and study volcanic activities, slow earthquakes, and tremors [1][2][3]. Moreover, strainmeters play an important role in detecting transient deformation associated with earthquakes and exploring whether atmospheric pressure variations in a typhoon are able to trigger fault slip [4][5][6][7]. Taiwan Island is located in a subtropical area with frequent typhoons characterized by significant air pressure drops of several tens of hectopascal (hPa) and heavy rainfalls resulting in geohazards [8,9]. Thus, typhoons represent one of the major hydrological factors triggering significant transient deformation in terms of atmospheric events in Taiwan [5,6].
Taiwan is located at the active orogenic belt between the Eurasian Plate and Philippine Sea Plate, with a convergence rate of 80 mm/yr [10][11][12][13][14]. The Philippine Sea Plate subducts northward beneath the Eurasian Plate in the northeast of Taiwan, forming the Ryukyu-arc Trench system. The continental crust of the Eurasian Plate subducts eastward beneath the oceanic crust of the Philippine Sea Plate in Southeastern Taiwan, forming the Luzon  The period of aseismic strain during the seismogenic process varies between one tenth of a second and several decades. Various detection techniques, such as seismographs, GPS, and InSAR (interferometric synthetic aperture radar), have been applied to observe the surface deformation in Taiwan [25][26][27][28][29][30][31] (deeply buried strainmeters are characterized by high precision and wide bandwidth; they directly record crustal deformation and efficiently reduce the effect of man-made noise above the ground [32,33]). Due to the high sensitivity of the borehole strainmeter, the data not only contains faint signals of seismic activity, but also the response to external environmental changes such as air pressure and rainfall will be found in the data. Figure 2 shows the observation data of the CINT station in 2011. The influence of rainfall can be detected by removing the responses of the air pressure and tide (arrows). The areal strain is compressive during periods with heavy rainfall. The earthquake event catalog shows there is no direct relation with these responses. The effect of rainfall on the strain has been frequently mentioned in groundwater observation studies [34,35]. It has been suggested that well-confined aquifers can be regarded as quasistrainmeters. In a number of studies, the crustal deformation could be calculated using the water level of the confined aquifer [36,37]. Based on poroelastic constitutive relations [38], the ratio between the water-level change of the confined aquifer and the crustal deformation could be calculated using the earth tide response. The motivation of this paper is that we want to characterize some abnormal strain regarding to the loading of rainfall observed in the borehole strainmeter. We try to understand their mechanisms, and seek numerical solutions to eliminate such anomalies, as the basis for further strainmeter data processing for the accumulation of tectonic strain in the fold-and-thrust belt of Taiwan.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 4 of 24 Figure 2. Responses of air pressure and rainfall observed on channel 3 (CH3) at station CINT. Light blue lines represent air pressure records; green lines are hourly rainfall; the light gray and purple bars indicate seismic events within 100 km and 50 km from CINT, respectively. ε is the residual strain removing atmosphere and tidal response, separately. After removing the response of air pressure and tides, the influence of rainfall becomes apparent (pointed by red arrows).

Strain Transformation and Instrument Coupling
The GTSM strainmeter can describe the strain state of an object under stress by calculating the data obtained from more than three transducers in different orientations [32], in which the extensile strain is set as a positive value and the compressive strain is set as a negative value in the original data. At a point in the body, let the Cartesian components Responses of air pressure and rainfall observed on channel 3 (CH3) at station CINT. Light blue lines represent air pressure records; green lines are hourly rainfall; the light gray and purple bars indicate seismic events within 100 km and 50 km from CINT, respectively. ε is the residual strain removing atmosphere and tidal response, separately. After removing the response of air pressure and tides, the influence of rainfall becomes apparent (pointed by red arrows).

Strain Transformation and Instrument Coupling
The GTSM strainmeter can describe the strain state of an object under stress by calculating the data obtained from more than three transducers in different orientations [32], in which the extensile strain is set as a positive value and the compressive strain is set as a negative value in the original data. At a point in the body, let the Cartesian components of horizontal displacement be u in the x direction (west-east) and v in the y direction (southnorth). In that case, the strain tensor can be expressed as three independent components, they are ε ee = du/dx, ε nn = dv/dy, and ε en = ε ne = 1/2(du/dx + dv/dy), respectively. Based on the definitions above, the engineering strains ε a , ε d , and ε s describing the strain state can be further deduced. The areal strain ε a is equal to ε ee + ε nn , which reflects the areal change; ε d is named the differential shear strain or γ 1 strain, equal to ε ee − ε nn , which are the west-east tensile behavior and south-north compressive behavior, respectively; ε s is defined as the simple shear strain or γ 2 , which is equal to 2ε en . The deformation behavior of a plane can be described using these three engineering strains [39,40].
Although the borehole strainmeter is stacked with multiple transducers in different orientations, compared with the scale of the installation depth, we still regard it as a plane strain. Thus, instrument strains recorded in different orientations can be expressed as engineering strains by the following relation: e i is the instrument strain in the axial orientation of ϕ i , which shows the angle counterclockwise from the x axis. However, due to such a small-scale geological inhomogeneity or installation situation, a cross-coupling terms is introduced into the expressions relating remote strains to instrument strains. For isotropic coupling, instrument and remote strain can be described as the following relation [41]: Although we generally do not consider the influences of vertical coupling, some cases were reported in some situations [40]. If vertical coupling is taken into account, the relation between the instrument and remote strain can be expressed as the following: where g i is the mechanical gain, ϕ i is for the angle counterclockwise from the x axis; C, D, and H represents axial, perpendicular, and shearing strain deformational parameters in the ith channel, respectively; Vε zz shows the effect of vertical strain where V is the vertical coupling parameter; and ε zz represents the vertical strain. The responded strain by vertical stress such as rainfall and barometric has the relation with the horizontal areal strain as the following equation described: In the equation above, ν is denoted as the Poisson ratio of the rock formation. Substituting Equation (4) into Equation (3) yields C is an apparent areal strain, which present the comprehensive result of areal and vertical strain, and is defined as the following: Appl. Sci. 2021, 11, 1301

of 21
As Equation (6) denotes, areal strain will be reduced by vertical coupling effects, apparent areal strain may vanish or become negative when V > (ν/(1 − ν))C in some location with a large vertical coupling effect [40]. If we assume the mechanical properties of rock are homogeneous and isotropic rock, H could be equal to D. Equation (6) could be reformed: Additionally, Equation (7) could be easily generalized as the following: Equation (8) describes the general situation of coupling. In most cases, the linear strain e i measured by the ith channel of strainmeter can be expressed with Equation (9): [A] 3×i is a coupling matrix that contains instrumental, angular, formation coupling and topographic effects comprehensively. The inverse matrix [A] −1 3×i is the calibration matrix, which can transfer to remote strain from observed data channels. In the different constrain, the method to obtain the coupling matrix will be different. Since the significant responses of rainfall, air pressure, and river loading contribute in can be concluded as some kind of vertical strain, we treat the data calibration as in the isotropic rather than vertical condition in this study. We had a short discussion of difference between isotropic and vertical coupling in Section 6.2.

Effects of Rainfall and Groundwater on the Strain
The data detected by strainmeters in Taiwan contain some response strains regarding natural environmental factors such as air pressure, tide, and especially hydrology. In a subtropical environment with an annual rainfall of 2500 mm on average such as in Taiwan, the hydrology generally has a strong effect on the observed strain. The effects of continuous heavy rainfall on the observed strain can be divided into two parts: (1) A large volume of rainwater falls in a short period of time, causing an increase of the external pressure on the rocks due to the extra loading of the rainwater. The rainwater that cannot infiltrate into the ground becomes land-surface runoff and rapidly merges into creeks. This will respond rapidly in hours since rainfall began, which is called the rainfall quick response (RQR); and (2) water slowly sinks into the strata to complement the groundwater level, which will last for several days to a few weeks due to the low permeability of the strata. This is called rainfall slow response (RSR) or base flow.
Strainmeter is designed to record the deformation behavior of the rock under stress. However, because the steel cylindrical strainmeter does not allow fluid flow into or out of the cylinder, the response of strainmeter is the same as the strain in the rock only when the rock remains undrained [42]. Checks the strain data during heavy rain, it shows a compressional strain within a short time after the loading of heavy rain, then gradually expands as groundwater pressure increasing with time ( Figure 3). This is different from the change of groundwater pressure under stress. When the rock is subjected to external stress in instant, the water in the rock cannot be drained immediately because of the low permeability of the rock, which can be considered as the undrained condition. In the undrained case, when the whole rock is compressed, the water pressure inside the rock will increase and the extra pore pressure will squeeze the strainmeter inwards. Both the instrument strain and rock strain show compression during this period. Next, the process of rainwater permeating into strata is continuous and slow, extra pore pressure can be balanced as water flowing. In this drained condition, the strainmeter showed expansion as groundwater pressure rose. However, these changes only reflect the balance of internal stress in the strata, there is no actual change in the far-field strain.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 7 of 24 Figure 3. The diagram shows the relation between the strainmeter strain and groundwater pressure induced by rainfall. Data from the DARB station located in Southern Taiwan (location shown in Figure 1).
Considering that the observed strain is closely related to the groundwater changes, a widely used recursive digital filter technique in hydrology was adopted in this study to separate quick and slow rainfall responses. The basic equation for recursive digital filter simply expressed as the following [43]: where Fi and Fi−1 are the filtered responses and Hi and Hi−1 are the strain at the ith and (i − 1)th sampling instants, respectively; is the filter parameter. Under the condition of 1800 sps (seconds per sample) preset in this study, we gave = 0.996 into the Equation (10) to create a high-pass filter that can effectively filter out the low-frequency signals with a period of around five days. The high-frequency signals can then be separated from the observed strain to identify the rainfall quick responses (RQRs).
We could obtain the strain corresponding to the rainfall load from the observed strain. Furthermore, we need to acquire the rainfall impact function Ri(t), which describes the whole responding process of each single unit of rainfall falling on the ground before flowing out of the effective area around the station. Rr(τ) is the rainfall record in half an hour in this study. Set τ as the unspecified time change. The RQR, R(t) is the convoluted result of the rainfall and impact function, which can be expressed with the following equation: Considering that the observed strain is closely related to the groundwater changes, a widely used recursive digital filter technique in hydrology was adopted in this study to separate quick and slow rainfall responses. The basic equation for recursive digital filter simply expressed as the following [43]: where F i and F i−1 are the filtered responses and H i and H i−1 are the strain at the ith and (i − 1)th sampling instants, respectively; λ is the filter parameter. Under the condition of 1800 sps (seconds per sample) preset in this study, we gave λ = 0.996 into the Equation (10) to create a high-pass filter that can effectively filter out the low-frequency signals with a period of around five days. The high-frequency signals can then be separated from the observed strain to identify the rainfall quick responses (RQRs).
We could obtain the strain corresponding to the rainfall load from the observed strain. Furthermore, we need to acquire the rainfall impact function R i (t), which describes the whole responding process of each single unit of rainfall falling on the ground before flowing out of the effective area around the station. R r (τ) is the rainfall record in half an hour in Appl. Sci. 2021, 11, 1301 7 of 21 this study. Set τ as the unspecified time change. The RQR, R(t) is the convoluted result of the rainfall and impact function, which can be expressed with the following equation: (11) In order to calculate the effect of rainfall on the strainmeter, we divided rainfall into many rain units. When this rain unit falls on the ground, and its weight will cause regarding deformation to the rock bed, although its weight is extremely small. As the rain unit flow out to the periphery of the observation station for a certain distance due to gravity, the effect of the rain unit will gradually decrease to be negligible. This distance is called the effective distance of the observation station. When the rain unit falls on the station (the force is the greatest now) to the distance from the station to the effective distance, a continuous strain-time path curve will be recorded by the strainmeter, which we call the rainfall impact function. By multiplying the rainfall impact function and continuous rainfall records in the frequency domain (i.e., convolution), we could obtain the rainfall curve for the entire time period. Strain records the response of the rock to stress. Therefore, rainfall impact function also reflects the transmission of stress. Once the rainfall falls, the maximum response is reached in a very short time. As the rainwater dissipates by time, the impact of the station will gradually decrease. The schematic diagram of the rainfall impact function is shown in Figure S2.

Data and Calculation
The observed strain change is often influenced by many environmental factors [32,[39][40][41]. In additional to dominant trend of borehole relaxation and grout curing, the strain changes are associated with tidal response, barometric pressure, rainfall, groundwater fluctuation, seismic deformation, and interseismic strain accumulation due to tectonic loading. To better understand the tectonic process from the observed strain, the non-tectonic factors such as barometric loading and rainfall response should be removed from the strain data. The borehole relaxation and grout curing can be removed by fitting the model of two exponential attenuation functions. Together with the removal of non-tectonic effects on the strain response, the intensity of tectonic signals will be largely improved.
In removing the tidal response, we used a software of SPOTL (Some Programs for Ocean-Tide Loading), see [44] to calculate theoretical tidal constituents with NAO.99B [45] as global, and TWTIDE08 [46] as the regional model. A program of BayTap-G [47] was used to analyze and extract the tidal component out of the observed data. Observed strain can be disturbed by severe hydrologic changes, such as typhoon and flood. Therefore, in order to reduce noise and improve the readability of signal, we mainly chose the data in the period of dry season (Oct.-Apr.) to calculate tidal constituents. In addition, averaging the calculation results of several years can be used as a reference basis for removing tidal response. The changes of air pressure also influence the observation in some linear relation. A period of 45 days was set as the sampling window; the barometric pressure and strain were respectively band-pass-filtered to exclude frequencies outside 3-5 days band [40], because in this time period, semi and diurnal signal can be ignored in a better signal/noise ratio and a more accurate response can be obtained. The barometric response strain, nstrain/hPa can be obtained by comparing both. The observed strain is not only affected by the surface water load but also by deep groundwater variations. Currently, the relation between the areal strain and groundwater pressure can be roughly correlated using co-sited water pore pressure transducer. However, the generated strain and the water pressure variation usually make it difficult to identify the tectonic activity-related strain, if some observation stations encounter the problems of water pressure gauge failures. The present work mainly focuses on the rainfall impact in a short period of time; the relation between the groundwater pressure variation and areal strain will be discussed in future publications.
The data from ten observation stations are selected for the calculations except SLIN in which data has been affected by borehole adjustment severely so that it is not capable to be taken involved in this study. The topography of the areas where the stations is located is so diverse that we can discuss the impact of rainfall on the strain observed in different environments. For reducing the cross-influences of various factors, seismic effects have been excluded from the priority data. The basic information about the stations is listed in Table 1. Rainfall events with a cumulative rainfall over 20 mm were selected because a recognizable rainfall response can hardly be observed from the data of events with little rainfall. Moreover, the comprehensive effect of multiple rainfall events may be observed frequently because a large number of long-lasting continuous rainfall events generally occur in the rainy season from May to November, demonstrating the relatively strong noise during this period. Therefore, independent heavy rainfall events in the dry season (December to April) are the priority to be selected as the study samples. Subsequent to the preliminary selection of rainfall events, relatively strong noise signals with a fairly high standard deviation would be removed before calculating the average value to reduce the impact of noise on the observed data. In this case, the stability of the computed results can be ensured.
Comprehensively, we briefly summarized the data processing flow to following steps: • First, convert the original data recorded by the borehole strainmeter into the actual deformation base on the prior calibration results. • Second, remove the responses of air pressure, tides, and borehole-curing to obtain the residual strain.
• Then, according to the rainfall records, select the strain data corresponding to the heavy rainfall time in the residual strain, and try to use a recursive digital filter to separate the quick and slow response caused by rainfall.

•
Calculate the deconvolution results of RQR and rainfall to obtain the rainfall impact function, which shows the strain-time curve of every falling rainwater unit before flowing away.

•
After an event finished, we went to check the next rainfall event, and repeated the steps described above.
We then averaged all the results to obtain a general function and apply it to all time of the station. Since the quick response caused by the rainfall load is the compressive strain, all the strain mentioned in the text below were denoted compressive as negative in the plots.

Results
The sampling period of the observation stations is shown in Table 1. The analysis of the RQR corresponding to the rainfall load was conducted on relatively heavy rainfall events (cumulative rainfall >20 mm) during the sampling period. The analysis results for the stations are illustrated in Figure 4, in the left panel for each BSM displays the rainfall impact functions obtained in different periods and the average function, and simulated RQR obtained by convolution of the average rainfall impact function and rainfall are shown on the right. The vertical axes in Figure 4 display responding strains and the horizontal show the time. According to the altitude and topography in the neighboring regions of BSMs, we roughly divided them into 2 categories, plain/low terrain and hill/mountain area, and checked the responses of rainfall in different topographic conditions.   The category of plain and low-terrain means those set in the regions of low altitude and small terrain change. JING is located at the edge of the Taipei Basin with an altitude of only 19 m. The neighboring region has a flat terrain with only one 10 m high hill approximately 100 m west of the station and without significantly large rivers passing within 500 m. Overall, there were no notable terrain effects on the borehole strainmeter observations. The results show that the rainfall impact function of JING was fairly stable, demonstrating no remarkable noise-induced residual influence on the data. The intensity of the hetero-phase noise could be effectively reduced by calculating the average. The areal strain at JING rapidly reached a peak of 8.8 nanostrain/hPa a half hour after the rain started and then descended over time. It took about 26 h to decrease the strength of responding strain to 10% of the peak value. Another low-terrain BSM was SANS, which is located in a flat terrain of the river bank on the edge of the Taipei Basin. A creek runs northward in the east of SANS at a distance around 200 m away. The rainfall response of SANS was similar to JING. After the onset of rainfall, the response reached the maximum quickly and then reduced gradually over time. However, the rainfall decay time in SANS persisted longer than in JING, it took more than 95 h after the onset of rainfall.
Constrained by the limitation of lithology and topography for installation, most of BSMs installed by CGS were located in hilly areas in the Western Foothill of Taiwan ( Figure 1). Apart from the topographic relief, rivers were close to the observation stations, there were 6 BSMs siding by a notable river crossing in the vicinity of the hill or mountain area.
TAIS is a typical BSM installed in the mountainous area. It is located at a slope with an altitude of around 800 m with relatively high terrain changes. According to the results, it took 0.5-1 h to reach the maximum rainfall response and persisted to at least 178 h long after raining began. The rock formation where the borehole strainmeter installed was composed by alternations of sandstone and mudstone, Late Miocene, and was overburdened by a thick colluvium and weathering debris. Water level was recorded around 50 m beneath the ground surface. In the surrounding area of TAIS, some gullies were observed, but there was no obvious water system developed, which was not conductive to concentrate runoff. A similar topographic condition was observed in the vicinity of the TSUN station. This might be the reason why the rainfall decay time of these two stations persisted longer than the others because the rainwater was slowed down by the rugged terrain.
Nevertheless, not all BSMs installed in the hill or mountain area show a long persistence in response, e.g., BMMT. BMMT is located on a terrace top in the hilly region with a river 250 m west of the station; the altitude is approximately 195 m. The responding areal strain of BMMT reached a peak of 5.23 nanostrain/hPa, and gradually declined over time in half hour after the rain started. It is notable that the rainfall decay time of BMMT did not last as long as the BSMs in hilly area such as TAIS and TSUN. Since the network of streams will affect the disputation of surface runoff [48,49], the river meandering around the terrace where BMMT sit might be one of the important factors for dispersing over-land water.
In addition, we found some significant features in DARB. Taiwan is located on the plate boundary and experiences a dramatic rise of strata under the mutual compression of the plates. Within the distance of about 100 km, the terrain has risen to thousands of meters high from the west coastal plain to the central ridge. It is therefore difficult to avoid the installation of borehole strainmeters on high mountains. DARB is located on a terrace of the upstream of the Zengwun River in Southern Taiwan, with an altitude of 953 m and surrounded by a large number of mountains >1000 m. There are some rivers and creeks around the station. The observed strain data of DARB were highly sensitive to rainfall, and the strain values and variations differed from those of other stations. Thirty rainfall events from 2007 to 2012 were selected for rainfall impact function calculations ( Figure 5). The results were concentrated in two ranges, relatively small responding strains below 2.0 nanostrain/hPa and higher strains above 2.0 nanostrain/hPa with a notable delay. Based on the strain distributions, we divided the data into two groups for the discussion: group I, smaller responses (twenty rainfall events); and group II, larger responses (ten rainfall events). The computed average rainfall impact functions for the two groups were plotted in the same coordinate system to distinguish the potential differences of both. Apart from the maximum responses, the duration values were also different apparently. The areal strain of group I rapidly reached the maximum responding strain of 1.53 nanostrain/hPa half an hour after the rainfall started and then declined over time. In contrast, the areal strain of group II reached a peak of 5.33 nanostrain/hPa 25 h after the rainfall starts. Both delay time and strain values of group II were higher than those of group I. They show a similar strain variation in the first five hours and rainfall decay time of the responding strain, persisting for 4.5 days, which may be associated with the area of catchment.
around the station. The observed strain data of DARB were highly sensitive to rainfall, and the strain values and variations differed from those of other stations. Thirty rainfall events from 2007 to 2012 were selected for rainfall impact function calculations ( Figure 5). The results were concentrated in two ranges, relatively small responding strains below 2.0 nanostrain/hPa and higher strains above 2.0 nanostrain/hPa with a notable delay. Based on the strain distributions, we divided the data into two groups for the discussion: group I, smaller responses (twenty rainfall events); and group II, larger responses (ten rainfall events). The computed average rainfall impact functions for the two groups were plotted in the same coordinate system to distinguish the potential differences of both. Apart from the maximum responses, the duration values were also different apparently. The areal strain of group I rapidly reached the maximum responding strain of 1.53 nanostrain/hPa half an hour after the rainfall started and then declined over time. In contrast, the areal strain of group II reached a peak of 5.33 nanostrain/hPa 25 h after the rainfall starts. Both delay time and strain values of group II were higher than those of group I. They show a similar strain variation in the first five hours and rainfall decay time of the responding strain, persisting for 4.5 days, which may be associated with the area of catchment. Figure 5. Rainfall impact functions of the DARB station are divided into two groups based on the rainfall responses. Basemap shows the proportion of discrete data of Group II within certain timespan. Group I includes the runoff and riverhydrological impact caused by rainfall (brown circles) and group II comprises the impact of debris avalanches on the observation stations caused by heavy rainfall (blue squares). The data of group II were stacked; the unit amounts were counted and shown as the results on the base diagram in which the darker blue area contains 65% and the lighter one contains more than 90% discrete data. Basemap shows the proportion of discrete data of Group II within certain timespan. Group I includes the runoff and river-hydrological impact caused by rainfall (brown circles) and group II comprises the impact of debris avalanches on the observation stations caused by heavy rainfall (blue squares). The data of group II were stacked; the unit amounts were counted and shown as the results on the base diagram in which the darker blue area contains 65% and the lighter one contains more than 90% discrete data.
In the previous session, we used a technique of deconvolution to estimate the individual rainfall impact function of BSMs. With stacking up the results of rain impact function obtained in a different time period, we could extract the average function for all events. In the case, it will strengthen the common and weaken the uncommon parts. However, the average functions often show discretely and noisily. Thus, we tried to find a function that could describe the behaviors of the most of average function in the continuous curve. In some studies, the runoff effect induced by rainfall could be fitted by an exponential curve under the limit of the power law [50][51][52][53]. According to the results, the response will reach the maximum in 1-2 h, then show an exponential decay with time. It can be described as the following equation: where h is the function of rainfall response changing with time. The coefficients α and β are constant, and the lag time k i at which has a maximum value of h(k) are needed to be determined. Due to the sampling frequency, it is impossible to obtain the exact lag time under half an hour. Since the response of all stations quickly reached the peak after the rain, all k i were set to be one. We calculated and concluded the results of all BSMs in Table 2 and all fitted curves are plotted in Figure 6. More, the response of RNT and RST was smaller than the others, and it might be because a different BSM model, GTSM100, was used here. It is the earlier developed version than GTSM21, and is reported to have low sensitivity to the vertical strain in earlier studies [40].  Base on Equation (12), the constant α means the maximum response at the time of ki, and β describe the behavior of strain decay by time. According to the results, all maximum responses were under ten nanostrains per hectopascal and were different individually. Base on Equation (12), the constant α means the maximum response at the time of k i , and β describe the behavior of strain decay by time. According to the results, all maximum responses were under ten nanostrains per hectopascal and were different individually. There was no significant relation between the maximum response and elasticity, the reason causing such a difference was not clear. Besides, the decay constant, β describes the time that a rainfall impact persisted. The quick response of rainfall mainly reflected the additional strain caused by the rainwater, which had not infiltrated into the groundwater after rainfall. The influence of runoff was declining gradually as the rainwater flowed out of the adjacent region of the station. One of the major factors for influence of runoff was terrain roughness, which involves the process of generation, convergence, channelization, and concentration of runoff [54].
The terrain roughness is generally defined as the ratio of the surface area of a surface unit to its projected area on a horizontal plane [55]. It is an important quantitative indicator to measure the degree of surface erosion. In some studies, terrain roughness has a certain impact on the collection and dissipation of surface runoff [56,57]. We tried to calculate the overall terrain roughness in the vicinity of each BSM to understand its relation with rainfall decay time. The digital terrain model (DTM) of 6 m resolution was taken to calculate the terrain roughness in the vicinity of 500 m radius of the BSMs, and the results are shown in Table 2. From the results, there seemed to exist a linear relation between the terrain roughness and rainfall decay time. However, we also found some BSMs in high roughness areas showing a relative short decay time. There were six BSMs (including BMMT, PFMT, CINT, DARB, RST, and RNT) located near rivers within a distance of 500 m in the mountain area, and rainwater could easily flow into the surrounding rivers. Thus, the quick dissipation of rainwater will result in the shorter decay time of the rainfall response (see Table 2 and Figure 6). For example, PFMT is located at a smooth slope of high altitude among mountains, which is conductive to concentrate the rainwater eastward to an existing river beside. Considering low terrain roughness and conduction of existing river, the rainfall decay time in PFMT took only 33 h to down to 10% of the initial response. Similar phenomena can also be observed in BMMT, RST, and RNT, which are expected having long-persisting decay time. On the other hands, the influence of rivers on runoff collection in plain/low areas seems to be smaller than that in hilly/mountain areas. For SANS, the terrain roughness is not high, and there are rivers passing nearby. The rainfall decay time of this station is expected to be short, but in fact its decay time was relatively long. The reduction of river pooling seemed to be not effective on the plains, or it might be affected by other factors. In addition, JING, which had the shortest decay time, was located on the plain and had a low terrain roughness. Sitting at the edge of the urban area and better artificial drainage facilities might be beneficial to shorten the decay time.
There are many factors affecting surface runoff flow rate, including slope, surface vegetation, biological activity, soil distribution, particle composition, soil saturation, and topographic fluctuations [54]. More studies are required to understand the mechanism of generation and dispersion of runoff.

External Factors Influencing the Rainfall-Induced Response at DARB
The previous sections showed that the results for DARB were not only rainfall-induced response, but also impacted by other factors, which might be attributed to its specific geological location. The observed data was therefore not only influenced by the river water volume, some other external forces such as debris avalanches might take an important role in the abnormal strain response. The results for DARB show that the maximum responses and the delay time of group II were larger than those of group I indicating that some additional forces were acting. The DARB station is located in Alishan Village, where it is an area with a high potential for debris flow disasters after heavy rain. The debris likely falls into rivers, adding to the bed load. This might be the reason for the reported phenomenon. The comparison between the cumulative rainfalls of all selected rainfall events with the data of groups I and II (Figure 7) indicates a threshold of 250 mm, which is essentially equal to the rainfall triggering index (RTI) of Alishan Village proclaimed by the Soil and Water Conservation Bureau of Taiwan (SWCB). Figure 7 shows that there is an extremely high possibility of debris falling into rivers together with the abundant rainfall when the cumulative rainfall exceeded 250 mm, forming debris flow in the river courses and generally impacting the observations. Although debris avalanches are closely related to rainfall, it is fairly difficult to accurately predict the occurrence time, location, and volume of debris, which substantially affects the delay time and maximum response of group II. Nevertheless, there are several strain curves at different times roughly have similar delay times and peaks, which might imply the repeated landslide happened at the same place.

The Influence of Vertical Coupling in Rainfall and Barometric Responses
In many studies, it has been mentioned that BSM data is calibrated by theoretical tides [5,32,[39][40][41][58][59][60]. Hart et al. proposed a calibration method in isotropic condition [41], Roeloffs [40] found some BSMs data set by the Plate Boundary Observatory (PBO) have strong vertical coupling, which will strengthen the response for rainfall or snow the surface [40]; Hodgkinson et al. believes that the orientation of gauges often erroneous due to the magnetic properties of rock or malfunction of the instrument compass, the proper adjustment of the initial orientation is suggested [39]. In our study, we can find some strong influence of vertical strain in the barometric response. We can effectively diminish the vertical strain by using Hodkingson's method, but the major purpose of this paper is focused on the dealing with those vertical strain caused by the surface loading such as rainfall, air pressure, and river loads. Therefore, we used the isotropic coupling method to correct the observations and highlight the vertical parts of the data. We compare two calibration methods of different considerations to see its difference.
Comparing the results of two different ways of calibration, the two tend to be con- Figure 7. Comparison of the cumulative rainfall of each rainfall event with the data of the two groups at the DARB station. The threshold is 250 mm, which shows that landslides are likely to happen when the cumulative rainfall is larger than 250 mm, increasing the sediment volume in the rivers. The threshold is in good agreement with the rainfall triggering index (RTI) (250 mm) of this area published by the Soil and Water Conservation Bureau of Taiwan (SWCB) of Taiwan.

The Influence of Vertical Coupling in Rainfall and Barometric Responses
In many studies, it has been mentioned that BSM data is calibrated by theoretical tides [5,32,[39][40][41][58][59][60]. Hart et al. proposed a calibration method in isotropic condition [41], Roeloffs [40] found some BSMs data set by the Plate Boundary Observatory (PBO) have strong vertical coupling, which will strengthen the response for rainfall or snow the surface [40]; Hodgkinson et al. believes that the orientation of gauges often erroneous due to the magnetic properties of rock or malfunction of the instrument compass, the proper adjustment of the initial orientation is suggested [39]. In our study, we can find some strong influence of vertical strain in the barometric response. We can effectively diminish the vertical strain by using Hodkingson's method, but the major purpose of this paper is focused on the dealing with those vertical strain caused by the surface loading such as rainfall, air pressure, and river loads. Therefore, we used the isotropic coupling method to correct the observations and highlight the vertical parts of the data. We compare two calibration methods of different considerations to see its difference.
Comparing the results of two different ways of calibration, the two tend to be consistent in the long-term but there may be some differences in the short-term, these differences may show the magnitude of the vertical strain. Taking the results of real strain in JING as an example, Figure 8 shows the comparison of corrected results and rainfall. The strain drop responding to rainfall events are still spotted on the results calibrated by the method of isotropic coupling. According to the results, the response of barometric change was around 10 nanostrain drops compared to a 1 hPa rise in air pressure. Nevertheless, it was about 7 nanostrain drops in the same situation in the results calibrated by the method of vertical coupling, decreasing the ratio about 30%. In addition, some rainfall response was weakened or even difficult to be observed after vertical coupling correction, and this will be not conductive the purpose of discussion of rainfall or river loads on the strain observation in this paper, so the method of isotropic coupling was preferred. According to the results, this response caused by local and external loads just can be observed in the sallow depth underground, it did not affect the strain state in the deeper strata in ten or even tens of kilometers, such as non-tectonic signals can be effectively reduced or even eliminated by appropriate correction methods, and the identifiability of the tectonic activity will be increased. was weakened or even difficult to be observed after vertical coupling correction, and this will be not conductive the purpose of discussion of rainfall or river loads on the strain observation in this paper, so the method of isotropic coupling was preferred. According to the results, this response caused by local and external loads just can be observed in the sallow depth underground, it did not affect the strain state in the deeper strata in ten or even tens of kilometers, such as non-tectonic signals can be effectively reduced or even eliminated by appropriate correction methods, and the identifiability of the tectonic activity will be increased. In some cases, the rainfall responses clearly spotted in isotropic-coupling data show weaker or even difficult to be identified after vertical coupling calibration.

Comparison of the Effect of Barometric Pressure and Rainfall
External stress will lead to the strain of strata. The barometric pressure and rainfall both are external forces but have different effects. The effect of barometric pressure is comprehensive and does not vary with the distance. In contrast, the effect of rainfall will gradually decline with the land-surface runoff flowing out of the range. However, both vertical coupling method, respectively. The dark colored line means the result of removing tidal response, while lighter color means no deal with tides. In some cases, the rainfall responses clearly spotted in isotropic-coupling data show weaker or even difficult to be identified after vertical coupling calibration.

Comparison of the Effect of Barometric Pressure and Rainfall
External stress will lead to the strain of strata. The barometric pressure and rainfall both are external forces but have different effects. The effect of barometric pressure is comprehensive and does not vary with the distance. In contrast, the effect of rainfall will gradually decline with the land-surface runoff flowing out of the range. However, both are induced by external stress, which should theoretically have an identical response under the same stress. By comparing the calculation of the areal strain of the barometric response and peak strain of rainfall impact (Table 2), the results demonstrated a linear correlation with a slope of 1 between the strains of the two factors (Figure 9), which was in good agreement with our assumption (Table 3). . Comparison between the barometric responding strain and maximum rainfall responding strain, which are approximately in phase.

Impact of Additional River Load on the Observed Strain
If an observation station is too close to a river, an additional strain will be generated under the impact of temporary water or debris weight, which is similar to a strip load placing by the station. Those loadings will exert an additional stress on the underlying strata. The strain detected by the strainmeter will decrease with increasing distance to a river. With the Young modulus (E) obtained from the uniaxial test, the additional strain

Impact of Additional River Load on the Observed Strain
If an observation station is too close to a river, an additional strain will be generated under the impact of temporary water or debris weight, which is similar to a strip load placing by the station. Those loadings will exert an additional stress on the underlying strata. The strain detected by the strainmeter will decrease with increasing distance to a river. With the Young modulus (E) obtained from the uniaxial test, the additional strain generated by the accumulated water or debris in the river bank after rainfall can be derived by calculating the additional stress on the strata caused by the vertical strip load [61]. It is known that the elastic modulus of the strata with DARB strainmeters was 4.25 GPa, the Poisson's ratio was 0.25 (assumed), the width of the river course was approximately 50 m ( Figure S3), and the debris density was given as 2680 kg/m 3 . Therefore, every 1 m think of debris accumulated in the river will generate a 0.1 microstrain response contributing to strainmeter (Figure 10). The observed data show that the severe Typhoon Fanapi that passed through central Taiwan in September 2010 caused 0.1 microstrain of compression in observation. Based on the computation results, it can be estimated that at least 1 m thick debris were produced in the Zengwun River bank near the DARB station during this typhoon. If the debris fell into the river along with the large volume of rainwater from the upstream (with given density of the earth flow 1200 kg/m 3 ), a 3-4 m thick earth flow would roll in the bank, which is equivalent to an earth volume >30,000 m 3 with a mass around 100 thousands ton. This would have significantly impacted the strainmeters near the river course. The area where DARB is located is a highly sensitive landslide area. Large-scale landslides may be triggered by heavy rain and typhoon every year. Since the DARB area is located at the meandering of the river, the debris from upstream will temporarily be deposited in the river bed. After many rainfalls, the deposited debris and rock will gradually be removed and carried downstream. According to the literature reported, the output of sand from Zengwun River where DARB is located upstream could reach 25 million tons every year [21]. We could find some terrace higher than 6 m besides the river, which indicates the essential deposited level of debris after earthflow ( Figure S4).
Appl. Sci. 2021, 11, x FOR PEER REVIEW 20 of 25 generated by the accumulated water or debris in the river bank after rainfall can be derived by calculating the additional stress on the strata caused by the vertical strip load [61]. It is known that the elastic modulus of the strata with DARB strainmeters was 4.25 GPa, the Poisson's ratio was 0.25 (assumed), the width of the river course was approximately 50 m ( Figure S3), and the debris density was given as 2680 kg/m 3 . Therefore, every 1 m think of debris accumulated in the river will generate a 0.1 microstrain response contributing to strainmeter ( Figure 10). The observed data show that the severe Typhoon Fanapi that passed through central Taiwan in September 2010 caused 0.1 microstrain of compression in observation. Based on the computation results, it can be estimated that at least 1 m thick debris were produced in the Zengwun River bank near the DARB station during this typhoon. If the debris fell into the river along with the large volume of rainwater from the upstream (with given density of the earth flow 1200 kg/m 3 ), a 3-4 m thick earth flow would roll in the bank, which is equivalent to an earth volume >30,000 m 3 with a mass around 100 thousands ton. This would have significantly impacted the strainmeters near the river course. The area where DARB is located is a highly sensitive landslide area. Large-scale landslides may be triggered by heavy rain and typhoon every year.
Since the DARB area is located at the meandering of the river, the debris from upstream will temporarily be deposited in the river bed. After many rainfalls, the deposited debris and rock will gradually be removed and carried downstream. According to the literature reported, the output of sand from Zengwun River where DARB is located upstream could reach 25 million tons every year [21]. We could find some terrace higher than 6 m besides the river, which indicates the essential deposited level of debris after earthflow ( Figure  S4).

Conclusions
A conclusion can be drawn based on the comprehensive results of the observed data from ten stations: regardless of the effect of the base flow on the areal strain, the rainfall will have a different impact on the RQR based on the geological and topographic differences of the surrounding areas of the stations. The rainwater that is not able to infiltrate into the ground immediately after falling merges with the land-surface runoff and induce the quick response of strainmeter. It demonstrates the maximum response in a short period of time due to the impact of the cumulative weight of the water and then exponentially decreased as rainwater flowed out of the range. The fluctuations of topography will affect the rainfall decay time, whereas the river beside might help accelerating runoff dispersion. The larger cumulative rainfall will have a stronger effect on rock flow-sensitive areas. When the cumulative rainfall reached a specific threshold in such an area, the whole area likely experienced a rainfall-induced landslide or debris flow. The responding strain curve potentially shows abnormal under the domination of landslide events.
In the long term, the strain trend will be affected by groundwater flow, but this process will range from 1 to 7 days, which was longer than surface runoff (Table 2). At present, we only scoped on the impact of surface overload on deformation observations in this study, and did not yet consider the effects of groundwater. To understand the whole proceeding of rainwater's influence on and under the ground, it is required to involve in the studies of underground hydrology and hydrodynamics. That will be important for our future research.
On the other hands, the observed results implied that the noise level will be increasing if a borehole strainmeter station settled by the river due to river hydrology and landslides. More effort will be required to reduce the noise level to acquire accurate observation results. Therefore, building a station in such an environment is not advisable. However, stations by the river courses are sometimes unavoidable for crustal strain observations due to the elevated terrain at the fold-and-thrust belt in ongoing mountain building of Taiwan orogenic belt. For the stations vulnerable to the impact of river flow variation, more careful considerations are required with respect to data processing and interpretation.
Supplementary Materials: The following are available online at https://www.mdpi.com/2076-341 7/11/3/1301/s1. Figure S1: The diagrams show the relief maps of the vicinity of borehole strainmeter stations. Figure S2: diagram of rainfall impact function. Figure S3: The diagram shows the relation of a specific underground point, P and a strip load with width of B past aside on the ground. Figure S4: The photo was taken on the riverside by the DARB.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.