Self‐Potential Ambient Noise and Spectral Relationship With Urbanization, Seismicity, and Strain Rate Revealed via the Taiwan Geoelectric Monitoring Network

Geoelectric self‐potential (SP) signals are sensitive to natural and anthropogenic factors. The SP spectral characteristics under the different factors in Taiwan were investigated, and the SP spectral scalings were correlated with urbanization level, seismicity, and crustal deformation. The ambient SP noise models were first established by estimating the probability density functions of the spectrograms at each frequency. The effects of the natural and anthropogenic factors on the SP signals are understood by comparing the SP noise models under various conditions, such as precipitation, urbanization, and electric trains. Results show that the SP signals in areas of high industrialization and human activity and areas close to train stations behave as white noises and exhibit a distinct spectral ripple at frequencies around 1 Hz. On the other hand, the SP spectral power law parameters, Gutenberg‐Richter b values, and dilation strain rates were estimated by using the SP, earthquake catalog, and GPS data, respectively, during 2012–2017. By investigating the correlations of the SP spectral parameters with the Gutenberg‐Richter b value, dilation strain rates, and urbanization level, the SP optimal frequency band is found between 0.006 and 1 Hz due to the high correlation between the SP and seismicity data and between the SP and dilation data and the low correlation between the SP and urbanization data. Hence, this study may help the filtering and screening of the SP data and facilitate the understanding of the mechano‐electric behavior in the crust.

The SP signals clearly exhibit daily variations, which are affected by tidal forces and human activities (Huang & Liu, 2006;Szarka, 1988;Uyeda et al., 2002). Hence, it is important to first understand the daily background noises, called ambient noises hereafter. However, it is extremely challenging to quantify the degree of human activities in certain areas, particularly for the degree of electricity usage, which is directly associated with the SP signals. Moreover, the SP signals are reported to behave as 1/f β noises (Balasco et al., 2002;Ramírez-Rojas et al., 2011). The Earth's crust is regarded as a self-organized critical system (Rundle et al., 2000, and references cited therein), exhibiting the fractal behaviors in mechanical and electromagnetic ways (Chen et al., 2006;Potirakis et al., 2017;Sornette & Sornette, 1990). Furthermore, crustal deformations are also suggested to affect the fractal behaviors (Öncel & Wilson, 2004). Hence, it is attracting to study the relationships among the strain rates and scaling factors of the SP and seismicity data.
In this paper, we took Taiwan as an example and investigated correlations among the data sets of the SP, seismicity, strain rates, and urbanization. At present, the anthropogenic factor was provided with a detailed spatial distribution in Taiwan only for the urbanization data (Huang et al., 2018), which may relate to human activities regarding electricity usage. Taiwan is located at the collision boundary between the Eurasian Plate and the Philippine Sea Plate at a convergence rate of 82 mm/year (Yu et al., 1997); hence, Taiwan has complicated geological features, active orogeny, and frequent earthquakes and provides multiple data sets concerning the SP, earthquake catalog, GPS (Global Positioning System), and urbanization data.
Since 2012, the Geoelectric Monitoring System (GEMS) has been established in Taiwan. It includes 20 SP stations about 50 km apart ( Figure 1). One of the main purposes of the GEMS network is studying correlations between SP signals and crustal activities (Chen et al., 2017;Chen & Chen, 2016;Jiang et al., 2019;Telesca et al., 2014). For instance, the SP anomalies were found to be associated with M L ≥ 5 earthquakes by using skewness and kurtosis and were attempted to create an earthquake forecasting algorithm (Chen et al., 2017;Chen & Chen, 2016). Moreover, associations between the time dynamics of the SP signals and crustal deformation processes were investigated by using the Fisher-Shannon method (Telesca et al., 2014). Obvious relations were found between the information properties of the SP signals and the differential strain intensity along the Taiwan orogen. However, due to a lack of ambient noise studies, it still remains contentious whether the GEMS network in Taiwan is producing meaningful SP data or this database is sufficient to study anomalous internally sourced fields.
In this study, we have characterized the SP 1/f β noises by investigating the continuous data for 6 years from the GEMS network. There are two main objectives: (1) to establish the ambient SP noise models for each station by estimating the probability density functions of the spectrograms at each frequency (Berger et al., 2004;Marzorati & Bindi, 2006;McNamara & Buland, 2004) and (2) to investigate the correlations of the SP spectral power law models with the urbanization level, seismicity, and dilation strain rate. We then discussed the behaviors of the SP noise models under several conditions, such as urbanization, electric trains, and precipitation. Also, we investigated the frequency dependence of the correlations of the SP spectral power law parameters with the urbanization level, seismicity, and dilation strain rate. To our knowledge, despite a few studies concerning long-term electromagnetic data (Han et al., 2016;Huang, 2011;Xu et al., 2013), no research has focused on the ambient noises and spectral scaling correlations of SP data from a network with densely distributed, similarly equipped, and long-term recording stations. Consequently, this study may contribute to enrich the understanding of mechano-electric features in a complex environment like Taiwan.

GEMS and SP Data
Since 2012, the GEMS network has been progressively installed in Taiwan and has been working by 20 SP stations uniformly distributed with a spacing of about 50 km (Figure 1). The SP data at each station are continuously recorded at a sampling rate of 15 samples per second and with a GPS time correction. The SP signals are the passive potential differences between two dipoles oriented in the north-south (NS) and east-west (EW) directions. The nonpolarized electrodes (Pb-PbCl 2 ) of each dipole were buried at a depth of 1 m to avoid ground disturbances (e.g., surface temperature, precipitation, and human activity) and connected with the local telephone cable to each recording office. Table 1 lists the name, location, and start date of each station. The information on the GEMS network can also be referred to a study (Chen & Chen, 2016).
In this paper, we analyzed the SP data for all components and stations measured from the onset time until 31 December 2017. Based on the 6-year data, we developed robust noise models for each SP station. To improve our understanding of the geoelectric features, we analyzed the daytime (09:00-14:00 local time, denoted as DT) and the nighttime (00:00-05:00 local time, denoted as NT) data and compared them with the whole day (00:00-24:00 local time, denoted as WD) data. In contrast to the DT data, the effects of electric trains and most human activities can be well avoided for the NT data.

Urbanization Data
The SP signals are sensitive to electric trains, factories, power pipelines, and other human activities concerning electricity usage (Fraser-Smith & Coates, 1978;Harada et al., 2004;Ishikawa et al., 2007;Oettinger et al., 2001;Saito et al., 2011;Takahashi et al., 2007). However, it is exceptionally challenging to quantify those above mentioned anthropogenic factors in certain regions. In this study, we employed the urbanization level estimated by Huang et al. (2018), who provided detailed spatial information on the degree of urbanization in Taiwan.
Urbanization is the process of the population moving from rural areas to urban areas, which is mainly proportional to the degree of human activities. Huang et al. (2018) studied the degree of urbanization in Taiwan townships and used data for the 2010 population, housing census, and so on. Four indicators in their study were introduced including "number of residents," "the percentage of people working in secondary industries," "the percentage of people working in tertiary industries," and "population density." Finally, the urbanization level was divided into five levels. Level 1 is the lowest, which had the lowest population size and density but the highest primary industry percentage. Level 5 is the highest, which had the highest population size, population density, and the highest tertiary industry percentage. Detailed analyses and categories of the urbanization levels can be referred to the study (Huang et al., 2018). The spatial distribution of the urbanization levels is mapped in Figure 1. The urbanization level of each SP station is also summarized in Table 1. However, the urbanization level here is neither time dependent during 2012-2017 nor directly related to the degree of electricity usage. In this study, we made two assumptions for comparisons between the urbanization and SP data. First, according to the Taiwan government reports (https://www.stat.gov.tw/sitemap.asp), the social change and economic growth have been becoming slow for the past decade; hence, we considered that the urbanization data for 2010 can be compared with the average parameters estimated from the data *The first three columns are the name, location, and the start date of a station, respectively. The fourth column is the number of the segments for estimating the PSDs using the WD, DT, and NT data. The fifth column is the urbanization level (U) estimated by Huang et al. (2018). The sixth column is the slope of the Gutenberg-Richter distribution (b) estimated by the catalog selected during 2012-2017 and within a 50-km radius of each SP station. The seventh column is the dilation strain rate (D) in μstrain/year derived from the GPS data during 2012-2017. The b and D values are represented by mean within one standard deviation. Note that (1) the urbanization level at KUOL is set U = 3 due to its proximity between areas of U = 2 and U = 4; (2) similarly, we set U = 2 at PULI between areas of U = 1 and U = 3.
sets (e.g., SP, seismicity, and GPS) during 2012-2017. Second, we assumed that the urbanization level is proportional to the degree of the factors regarding electric trains, factories, and human activities regarding electricity usage. Therefore, the relationship between the electricity-related factors and SP data can be understood to some extent by studying the relationship between the urbanization and SP data.

Earthquake Data
From 1973 to 1992, the Taiwan Telemetered Seismographic Network (TTSN; Wang, 1989) was operated by the Institute of Earth Sciences, Academia Sinica. The TTSN consists of 24 stations, each equipped with a vertical high-gain and analog velocity seismometer. In 1992, the TTSN was merged into the Taiwan Seismic Network operated by the Central Weather Bureau (CWB) to form the CWB Seismic Network (CWBSN; Shin, 1992;Shin et al., 2013). The CWBSN is composed of 72 stations, each equipped with three-component digital velocity seismometers (Wu et al., 2008) and provided high-quality digital earthquake data that are recorded in both high-and low-gain formats. Since 2010, the CWBSN has again been upgraded: The seismograms are currently digitized at 24 bits, and the number of stations has increased with more than 100 now widely installed in Taiwan (Chang et al., 2012;Hsiao et al., 2011). The types of seismic instruments include short-period seismographs, accelerometers, and broadband seismometers. Therefore, the source parameters of an earthquake can be estimated more accurately and the completeness of magnitude (M c ) in this network can reach down to 1.5 (Chen et al., 2013).
That Earth's crust was observed to exhibit fractal behaviors in mechanical and electromagnetic ways (Chen et al., 2006;Potirakis et al., 2017). Hence, we paid attention to study the relationship between the scaling factors of the seismicity and SP data. In this study, we employed the Taiwanese earthquake catalog during 2012-2017 with M L ≥ 0 in longitudes 119.5-122.5°E and latitudes 21.5-25.5°N. We then estimated the parameters of the Gutenberg-Richter (GR) law (Gutenberg & Richter, 1944;Rabinovitch et al., 2001), such as b and M c values, by using the catalog selected within a 50-km radius for each SP station.

GPS and Data Acquisition and Processing
Crustal deformations affect the fractal behaviors in the crustal systems (Öncel & Wilson, 2004). Hence, it is necessary to study the influence of deformation rates on the SP and GR scaling factors. In this study, we estimated dilation strain rates by using GPS data. Until 2019, the GPSLAB (http://gps.earth.sinica.edu.tw) has been managed more than 500 GPS stations operating in Taiwan. The GPS data were employed to perform numerous studies regarding crustal deformations and plate motions in Taiwan (Takahashi et al., 2019;Yu et al., 1997Yu et al., , 1999Yu & Kuo, 2001).
Here, we analyzed the GPS data during 2012-2017 as the same observed period as the SP and earthquake data. The GPS data were downloaded from the GPSLAB and processed with the Jet Propulsion Laboratory's GIPSY-OASIS software (Bertiger et al., 2010;Zumberge et al., 1997). This software corrects the effects of the Jet Propulsion Laboratory's final orbit and clock information, the antenna phase center calibration, the first-(and higher-) order ionospheric bias, the tropospheric delays, the ocean loading, and so forth. These positions were computed under a new International GNSS Service (i.e., IGS08) reference frame (Rebischung et al., 2012). Then, we estimated velocities and velocity uncertainties from the GPS coordinate time series with the Median Interannual Difference Adjusted for Skewness algorithm (Blewitt et al., 2016). This approach computes the median of slopes between all data pairs separated by 1 year, which is a variation of the Theil-Sen nonparametric median trend estimator (Sen, 1968;Theil, 1992). The difference of pairs separated by 1 year mitigates the effect of annual periodicity rather than other transient signals such as different periodicities and coseismic and postseismic deformations. This technique has been shown to be more robust than the traditional least squares with seasonal fitting (Blewitt et al., 2016) against common error sources such as step discontinuities, outliers, and seasonality.

Methodology
The SP signals are affected by daily variations, such as solar-terrestrial interactions, electric trains, and human activities regarding electricity usage. In order to understand the influences of natural or artificial sources on the SP signals, it is important to investigate the ambient SP noises, which are the daily background noises. In this study, we first established the ambient noise models for each SP station by computing the probability density functions (PDFs) on a set of the power spectral densities (PSDs) of selected SP signals.
The methodology presented here was mainly referred to and modified from earlier studies (Berger et al., 2004;Marzorati & Bindi, 2006;McNamara & Buland, 2004), which applied the method for evaluating ambient earth noises. Hence, comparisons of the SP noise models can be carried out under several conditions, such as the urbanization level, precipitation, and electric trains (see section 4).
On the other hand, the Earth's crust is a complex self-organizing system, which is composed of multiple physical processes that interact between individuals to give rise to the overall state of the Earth as an entirety (Rundle et al., 2000;Turcotte, 1992). Hence, the fractal behaviors in seismicity and electromagnetic data can be observed (Chen et al., 2006;Potirakis et al., 2017). Also, the fractal behaviors are affected by crustal deformations (Öncel & Wilson, 2004). Hence, in the second part, we analyzed and investigated the features of the SP 1/f β noises, GR b values, and dilation strain rates. The processing steps and estimating methods are detailed as follows.

SP Data Processing
The SP signals each day for every component and station were selected into one analyzed segment considering three kinds of time spans: (i) the whole day (00:00-24:00 local time) as the WD data, (ii) the daytime (09:00-14:00 local time) as the DT data, and (iii) the nighttime (00:00-05:00 local time) as the NT data. The effects of electric trains and most human activities can be well avoided for the NT data as a quiet reference. Figure 2a illustrates the original SP signals on 1 January 2016 for both components at CHCH. The SP amplitudes usually vary from tens to hundreds of mV/km with some extremely spiky signals. To better analyze the background noises, we thus discarded the data points with values greater than the 99.5th percentile and smaller than the 0.5th percentile in each segment. After such first data screening, if the percentage of missing data is less than 30% of one analyzed segment (at which the PSDs can still be well estimated by the Lomb-Scargle method; Munteanu et al., 2016; see section 3.2), then the unevenly spaced time series are further processed. Next, we calculated the mean value and linear trend of each segment and subtracted them from that segment. The total number of the remaining segments for later estimating the PSDs using different time-spanned data at each station is listed in Table 1.

PSD
The PSDs of the segments for different time-spanned data were estimated by the Welch method (Seats et al., 2012;Welch, 1967) and the Lomb-Scargle method (Lomb, 1976;Munteanu et al., 2016;Scargle, 1982). In each segment, we calculated the PSD in a moving window length of 10,000 s by using the Lomb-Scargle method with a shifting length of 1,000 s. Then, we computed a modified periodogram by averaging the estimated Lomb-Scargle PSDs in the segment. This averaged periodogram with overlapped time windows tends to reduce the variance for a single periodogram. This moving window is windowed by the Tukey window with a cosine taper of 25% width to decrease sidelobe leakage effects and to enhance the signal-to-noise ratio (Bloomfield, 2004). The window length has 10,000 s so that the longest resolvable period is roughly 1,000 s (10 −3 Hz) following an empirical relationship that the reliable longest period is about one tenth of the calculating window (McNamara & Buland, 2004). The frequencies are evenly spaced in a logarithmic scale with a step of 0.02 between 10 −3 and 10 0.86 (≅ 7.24) Hz. We then applied one-third octave smoothing to the PSD of each segment for reducing measurement noise (Hatziantoniou & Mourjopoulos, 2000;O'Haver, 2016;Tylka et al., 2017). Let us take the SP time series shown in Figure 2a as an example. Figure 2b demonstrates the PSD using the WD data for the NS component estimated by the above mentioned method. We observed that the PSD between 0.001 and 0.1 Hz is approximately linear in a log-log scale. The linear fitting of the PSD was then estimated by using the least squares regression (see section 3.4). Furthermore, we estimated the daily PSDs using the different time-spanned data for all components and stations. Figure 3 illustrates the spectrogram at CHCH from its start date to 31 December 2017.

PDF
To investigate the ambient SP noises, we at each frequency calculated the PDF of the PSDs for the different time-spanned data. We grouped the PSD values into bins of 0.1 log 10 (mV 2 /km 2 /Hz) wide, ranging from −5 to 7, obtained histogram, and then divided by the total number of PSDs. Figures 4a and 4b demonstrate the PDFs for the WD data at frequencies of 0.003, 0.03, 0.3, and 3 Hz at CHCH. We observed that for both components the PDFs at lower frequencies, particularly when f = 0.003 Hz, are usually skewed to the high power densities. This implies that most of affecting factors on the SP signals seriously disturb the low-frequency signals. Figures 4c and 4d further demonstrate the probabilistic power spectral densities (PPSDs), which represent the PDFs for all frequencies in one plot. We also estimated and plotted the 5th percentile (P 5 ), mode (P M ), 95th percentile (P 95 ) at each frequency, which are referred to as the ambient SP noise models hereafter. In this case, the percentages around the P M (f) are greater than 12% for both components. Moreover, the distance of P 95 (f) − P 5 (f) decreases with the frequency. This suggests that the SP signals at the low frequencies (f < 0.1 Hz) are more sensitive than those at the high frequencies (f > 0.1 Hz).
Figures 5a and 5b show the PPSD images using the DT data at CHCH for the NS and EW components, respectively. In contrast, Figures 5c and 5d show the PPSD images using the NT data for the NS and EW components, respectively. The percentages on the P M (f) for the DT and NT data for both components are all greater than 12%. On the other hand, the distances of P 95 (f) − P 5 (f) with f < 0.1 Hz for the NT data are smaller than those for the DT data. The day-night difference in this case logically suggests less interference of artificial noises during the nighttime. The distances of P 95 (f) − P 5 (f) with f > 10 −1 Hz for the NT data for the EW component are wider than those for the NS component. This might be due to the induced telluric currents by the interference of the Earth's magnetic field or the telluric current circulation from land to sea (Simpson & Bahr, 2005;Zhang et al., 2019).

Scaling Exponent of SP 1/f β Noise
The time dynamics of observational signals in the crustal system may exhibit correlation structures and selfsimilarity in mechanical (Chen et al., 2006;Sarlis et al., 2015;Turcotte, 1999) The scaling exponent β is usually observed in the range 0 ≤ β ≤ 3 and is generally estimated as the slope of the line fitting in the least squares sense of equation (1) presented in log-log: where α is a constant.
As a rule of thumb, a candidate power law should meaningfully exhibit a linear relationship in a log-log scale over approximately 2 orders of magnitude in both axes (Stumpf & Porter, 2012). Hence, we selected a moving frequency window with 2 logarithmic orders with a shifting step of 0.2 logarithmic orders. For example, we fitted equation (2) to the PSDs between 10 −3 and 10 −1 Hz, between 10 −2.8 and 10 −0.8 Hz, and so on, until fitting between 10 −1.2 and 10 0.8 Hz. In this study, we attempted to understand whether there is a frequency-dependent relationship of the SP data with the urbanization data, GR b value, and dilation strain rate. Using the WD data by fitting equation (2) to the PSDs between 0.001 and 0.1 Hz, Figure 6 illustrates at CHCH the daily variations of the estimated β values and the coefficient of determination (R 2 ) that quantifies how good the regression line fits the data. We observed that most of the β values for both components are in the range of 1-2 with R 2 > 0.8.
Further, we calculated the mean values and standard deviations of the daily estimated β and α values with R 2 > 0.8 for each SP station during 2012-2017. Table 2 shows the mean values of β and α (denoted as β and α) and the standard deviations of them (denoted as σ β and σ α ) for all components and stations by using the different time-spanned data. In this case of fitting the PSDs between 0.001 and 0.1 Hz, Figures 7a and 7b also illustrate the spatial distributions of the β values using the NT data (denoted as β NT ) for the NS and EW components, respectively. We observed that for the different time-spanned data and both components, the β values for most of the SP stations range from 1 to 1.5 with σ β from 0.2 to 0.4. This indicates that the ambient SP noises are often characterized by the flicker-noise and Brownian-noise processes (Guzmán-Vargas et al., 2009;Ramirez-Rojas et al., 2004). On the other hand, the α values for most of the SP stations range from −2 to 0 with σ α from 0.5 to 1.5. However, the scientific community was less interested in the discussion of the α values. In this paper, we attempted to discuss the related issue regarding the α values in the SP noises with the urbanization level, GR b value, and dilation strain rate.

GR b Value
The frequency-magnitude distribution of seismicity is generally characterized by a power law behavior (Gutenberg & Richter, 1944), also called the GR law: where N(m ≥ M c ) is the number of earthquakes with magnitude (m) larger than or equal to the completeness of magnitude (M c ), b is a scaling parameter, and a is a constant. Earthquake catalogs are shown to be incomplete at small magnitudes due to a detection threshold of seismic instruments and networks (Rydelek & Sacks, 1989); hence, assessing the M c value is critical. In this study, we estimated the b and M c values by using the Clauset method (Clauset et al., 2009;Virkar & Clauset, 2014), which optimizes the lower data cutoff (i.e., completeness of magnitude in seismology) by minimizing the Kolmogorov-Smirnov distance between the empirical cumulative probability function obtained from the data and the analytical cumulative probability function obtained via the GR law parameters. The detailed application of the Clauset method can also be referred to a study (Li et al., 2015).
In this study, the b value was estimated by the earthquake catalog selected within a 50-km radius of each SP station, which is also the interstation distance. Table 1 shows the b value corresponding to each station. The errors in the b values were also calculated by using a bootstrap approach (Schorlemmer et al., 2003). In this  Table 2 for the (a) NS and (b) EW components. Maps of (c) the GR b values and (d) the dilation strain rates (D) listed in Table 1. approach, we randomly resampled a data set of the same number of events on the selected catalog for each station and used the Clauset method to estimate the b value of the resampled data set. We repeated this procedure 1,000 times and took the mean value and standard deviation to represent the confidence in b (denoted as b and σ b ). As Figure 7c shows the spatial distribution of the estimated b values, which are close to 1 for most areas. However, the middle and western sides of Taiwan are characterized by smaller b values (below 0.9) than the eastern side. Furthermore, the northern sides are featured by higher b values (above 1.1); particularly, the b value (b = 1.74 ± 0.29) is maximum at SHRL located on the Tatun Volcano Group, in agreement that volcanic areas are usually featured by higher b values (Pu et al., 2014;Wiemer & McNutt, 1997;Wyss et al., 1997).

Dilation Strain Rate
The dilation strain rate (D) is the first invariant of the strain rate tensor (_ ε ij ), which equals the trace of _ ε ij : Since the GPS does not provide reliable vertical velocities (Schmidt et al., 2003;Wing & Frank, 2011), we considered the areal dilation rate, which is the sum of the maximum and minimum principal strain rates (_ ε 1 and _ ε 2 , respectively): The positive and negative D values correspond, respectively, to extension or elongation and contraction or shortening. We derived strain rates from the GPS velocities based on the VISR (Velocity Interpolation to Strain Rates) method described by Shen et al. (2015). In this method, we meshed the Taiwan area and set a grid size of 0.05 × 0.05°. The local strain rate at each node of the grid was estimated from local geodetic measurements using Gaussian and Voronoi cell weightings to automatically reduce the contribution of far stations from the node within the entire network. A uniform strain rate field was assumed at each node of the grid, and a weighted least squares inversion was performed over strain component solutions and their covariance for the unknowns of translation, rotation, and strain. Next, we estimated the derived quantities of the strain rates such as the principal strain rates and dilation strain rates.
We then selected the nodes of the grid within a 20-km radius of each SP station and averaged the dilation strain rates of the selected nodes. Table 1 lists the mean values and standard deviations of the dilation strain rates (denoted as D and σ D ). Figure 7d shows the spatial distribution of the D values. The large positive D values (above 0.2 μstrain/year) feature around TOCH on account of the back-arc rifting or extension in the Okinawa Trough (Liu et al., 2016;Sibuet et al., 1987). On the other hand, the large negative D values (below −1.6 μstrain/year) appear at YULI and RUEY due to the strong convergence of the Luzon arc (Mouthereau et al., 2009;S. Yu et al., 1999).

Ambient SP Noise Models
In this section, we paid attention to the comparisons of the SP noise models (P 5 , P M , and P 95 ) under different conditions, such as the urbanization level, electric trains (associated with the day-night difference), and precipitation.
To understand the influence of the urbanization level on the SP data, we compared the SP noise models for the WD data at KAOH with those at HUAL, as shown in Figure 8. Station KAOH (U = 4) is situated in highly industrialized areas in Taiwan and close to the area of U = 5, while HUAL (U = 1) is located in desolate mountain areas as a silent spatial reference. We observed that the SP noise models for both components at KAOH are flatter than those at HUAL, which are close to white noises. Moreover, distinct ripples appear in the noise models at KAOH, particularly at frequencies of 0.01 and 1 Hz. These results indicate the serious interference by human activities at KAOH, the highly industrialized area.
Second, to understand the effect of the electric trains on the SP data, we compared the SP noise models for the DT data with those for the NT data at HERM and HUAL, as shown in Figure 9. Station HERM is close to 10.1029/2019JB018196

Journal of Geophysical Research: Solid Earth
a busy train station with a distance of 5 km, whereas there is no train effect at HUAL in remote mountains. Additionally, the DT data are spanned during the working hours of 09:00-14:00 local time; contrarily, the NT data are spanned during 00:00-05:00 local time, when the train effect can be totally avoided, as a quiet temporal reference. We observed that at HERM and for both components, the SP noise models when f > 10 −0.5 Hz during the nighttime are 0.5 log 10 (mV 2 /km 2 /Hz) smaller than those during the daytime. Exceptionally, the P 95 models around f ≅ 1 Hz feature a distinct ripple and have the same order for both DT and NT data. Moreover, the distances of P 95 (f) − P 5 (f) with f < 0.01 Hz during the nighttime are smaller than those during the daytime. As the quiet reference at HUAL, the SP noise models are overall similar for both the DT and NT data. Only the P 95 model for the NS component during the nighttime is smaller than that during the daytime.
Third, to understand the impact of precipitation on the SP data, we compared the SP noise models at the quiet station HUAL and using the WD data, which were reestimated by the proposed methodology under

Journal of Geophysical Research: Solid Earth
the conditions of precipitation equal to 0 mm/day, greater than 50 mm/day, and greater than 100 mm/day. Figure 10a shows the daily time series of precipitation at HUAL, which displays the seasonal variations that the precipitation amounts are highly focused in summer times. Figures 7b and 7c show the ambient noise models with precipitation equal to 0 mm/day and greater than 50 mm/day for the NS and EW components, respectively. Here the SP noise models with precipitation equal to 0 mm/day are served as a dry reference. Figures 7d and 7e also show the SP noise models with precipitation equal to 0 mm/day and greater than 100 mm/day for the NS and EW components, respectively. We observed that the P M models are almost the same for the three precipitation conditions. However, the P 5 model for the NS component with

Journal of Geophysical Research: Solid Earth
precipitation >100 mm/day has larger PSDs for all frequencies than that with precipitation >50 mm/day, which has also larger PSDs than that with precipitation of 0 mm/day. Moreover, the P 95 model for the EW component with precipitation >100 mm/day has larger PSD when f > 10 −0.7 Hz than that with precipitation >50 mm/day, which has also larger PSDs than that with precipitation of 0 mm/day. Overall, whether the precipitation is considered, the ambient SP noise models show little difference. The reason may be that long dipole lengths are used (1-3 km) in the GEMS network. While the precipitation can cause a prominent flow direction in small scales and the strong effect of the precipitation on the SP signals in short dipole lengths such as tens of meters, the GEMS data look immune to rainfalls.

Correlations Among SP Power Law Parameters, Urbanization Level, GR b Value, and Dilation Strain Rate
In this section, we focused on analyzing spatial correlations (symbolized as ρ hereafter) among the urbanization level (U), GR b value, dilation strain rate (D), and SP power law parameters (α and β) using the different time-spanned data and for both components. By the way, we simply denoted the correlation between U and b as ρ U¯b , for instance. Figure 11 illustrates the correlation and corresponding p value matrices when the α and β values were estimated in the frequency band of 0.001-0.1 Hz. The correlation matrices for the NS and EW components are imaged as Figures 11a and 11c, respectively; correspondingly, the p value matrices for the NS and EW components are imaged as Figures 11b and 11d, respectively. We observed that the absolutes of correlations with |ρ| > 0.1 have significance with p ≤ 0.1. The correlations of ρ αWD¯αDT , ρ αWD¯αNT , and ρ αDT¯αNT are highly correlated, and so are those of ρ βWD¯βDT , ρ βWD¯βNT , and ρ βDT¯βNT . This logically indicates that the SP power law models at each station have similar patterns whether they are estimated using the WD, DT, or NT data. No spatial correlation exists between U and D (ρ U¯D = 0.1), which is reasonable that the crustal deformations are mainly caused by plate motions and mantle convection (Pysklywec & Beaumont, 2004;S. Yu et al., 1999). However, a moderate correlation seems to exist between U and b ( ρ U¯b = 0.39; see section 5). Moreover, a moderate correlation might also exist between b and D (ρ b¯D = 0.41) in agreement with the finding that seismicity is more common at positive strain rate values (Hauksson, 2011;Nakamura & Kinjo, 2018).
We now concentrated on the correlations of the SP power law parameters with the urbanization level, GR b value, and dilation strain rate. In this case of α and β estimated between 0.001 and 0.1 Hz, we observed that the correlations between U and β for the different time-spanned data are featured by relatively high negative values (−0.65 to −0.47) for both components. This means that larger U values correspond to smaller β values, which is reminiscent of the case of KAOH (Figure 8). We considered that the SP signals are seriously dominated by white noises when an SP station is located in highly urbanized areas. However, in this case, the correlations of α and β with b and D almost do not exist (|ρ| ≤ 0.3).

Journal of Geophysical Research: Solid Earth
To further understand the influence of the urbanization level, GR b value, and dilation strain rate on the frequency-dependent SP power law parameters, we then analyzed the correlations of U, b, and D with α and β estimated in different frequency bands (see section 3.4). Figure 12 shows the correlations of U, b, and D with α estimated in different frequency bands for both components using the different time-spanned data. We observed that no matter what frequency bands are involved, most of them are not correlated (|ρ| ≤ 0.3). Exceptionally, the correlations between U and α DT at some frequency bands can be featured by large negative values (below −0.3). This again implies that high human activities during the daytime indeed seriously affect the SP signals. On the other hand, Figure 13 similarly shows the correlations of U, b, and D with β estimated in different frequency bands. We observed that at certain frequency bands, moderate correlations of U, b, and D with β can reach values of |ρ| > 0.3. This result indicates that the correlations are frequency dependent. Compared those curves, the optimal frequency bands (10 −2.2 to 1 Hz) can be found due to higher positive correlations of ρ b¯β and ρ D¯β and less correlation of ρ U¯β .

Discussion and Conclusions
The effects of anthropogenic factors (e.g., electric trains, factories, and power pipelines) are reported to disturb electromagnetic signals (Fraser-Smith & Coates, 1978;Ishikawa et al., 2007;Saito et al., 2011;Viljanen & Pirjola, 1994). However, it is difficult to quantify such anthropogenic factors in certain areas. Here we choose the detailed spatial distribution of the urbanization level in Taiwan estimated by Huang et al. (2018) on behalf of the anthropogenic factors. In this paper, we demonstrated the comparisons of the ambient SP noise models under different artificial factors, such as the urbanization level and electric trains. As Figure 8 shows, the comparison for the urbanization level suggests that the SP signals are characterized by white noise processes in areas of high industrialization or human activity. On the other hand, in Figure 9, the day-night difference, which mostly responds to the train effect, between the SP noise models indicates that the SP signals with f > 10 −0.5 (≅ 0.32) Hz and f < 0.01 Hz are mainly affected by the train noises during the daytime. Hence, installing one SP station in the above mentioned areas should be avoided. However, due to site limitations, the SP signals are sometimes inevitably disturbed by the anthropogenic factors. Logically, a further and critical step is to develop a noise reduction method or apply developed methods to distinguish noises from signals, such as remote reference (Larsen, 1997;Oettinger et al., 2001) and interstation transfer function (Harada et al., 2004;Takahashi et al., 2007). The ambient SP noise models proposed here play an important role to examine whether noises are successfully extracted from signals after using noise reduction methods. For example, we would expect that the SP noise models become flicker-noise or Brownian-noise patterns and the spectral ripples in the noise models disappear if one applied the interstation transfer function approach to the SP signals at KAOH. At present, this remains an open question and will be our future work.
The most important topic is to study the frequency dependence of the correlations of the SP power law parameters (α and β) with the urbanization level (U), GR b value, and dilation strain rate (D). As Figure 12 shows, the α values are not correlated with the U, b, and D values. However, the α represents the energy level of the analyzed signals; hence, we believed that the α values in the SP power law models should be related to the geological settings or hydrological conditions, which can affect the electromagnetic signals (Befus et al., 2014;Rizzo et al., 2004). On the other hand, the correlations of β with U, b, and D show the frequency dependence, and the optimal frequency band is found between 10 −2.2 (≅ 0.006) and 1 Hz, at which there are relatively high correlations between β and b and between β and D and relatively low correlation between β and U. The Earth's crust, as a complex dynamic system, exhibits the fractal behaviors in the mechanical and electromagnetic ways in the same time (Rabinovitch et al., 2001); hence, the fractals in the seismicity and SP data logically have a positive correlation, which reaches a maximum of 0.6 in this study. Moreover, the correlation between β and D reaches a maximum of 0.4. As Öncel and Wilson (2004) pointed out, the crustal deformations can affect the fractal behaviors in the crust. Therefore, there are positive correlations among β, b, and D in this study. However, it is still unsolved that the effects of the anthropogenic and crustal factors on the SP signals between 0.006 and 1 Hz are better performed.
In the investigation of the spatial correlation matrix, interestingly, the positive correlation between the urbanization level and GR b value seems to exist. The urbanization level is mainly related to the degree of population, which is related to the surface mass loading. Hence, the high urbanization level represents large surface mass loading. It has been reported that the surface mass loading may affect the seismicity, such as snow load (Heki, 2003) and reservoir mass (Mandal et al., 2005;Simpson, 1976). Additionally, the coverage and detection of the seismic network should be considered. The seismometers in Taiwan are densely distributed in cities (high urbanization level) rather than in mountain areas (low urbanization level; Shin et al., 2013); hence, more small earthquakes in highly urbanized areas may be recorded and the GR b values are then affected to increase. However, there is no rigorous analysis and inference in this paper; hence, it is still arguable that earthquakes may be triggered by the population mass loading or caused due to the detection capability of the seismic network.
In conclusion, we have investigated the ambient SP noise models under different conditions, such as the urbanization level, electric trains, and precipitation. Understanding the patterns of the ambient SP noises is helpful to maintain the original SP stations or deploy new stations. As a typical example, the SP signals at KAOH are atrocious. Second, we have investigated the spatial correlation matrix and its frequency dependence of the SP power law parameters with the urbanization level, GR b value, and dilation strain rate. The determination of the optimal frequency band allows us to filter and screen the SP signals and improve the quality of the SP analyses. Consequently, this paper may be conducive to understand the seismo-electric theory and mechano-electric features in the crust.