Fluctuation Characteristics of Water Level and Water Temperature of Huize Well Based on MF-DCCA

Subsurface fluids are closely related to the formation, occurrence, and development activities of earthquakes. Under the premise that the confined aquifers of groundwater are in an intact state, the subsurface fluids can flexibly and quickly respond to changes in stress and strain of the earth crust. After the Xingtai earthquake in 1966, China began to establish observation stations to monitor groundwater conditions, after decades of continuous upgrade and perfection, now a subsurface fluid observatory network that can monitor the well water level, temperature, and geochemistry has been constructed. The Earthquake Cases of China recorded that, at the time of the 2004 M8.7 Sumatra earthquake, the water level of 78 wells and the water temperature of 59 wells in China showed co-seismic changes [1]. After researching the co-seismic response mechanism of water level and water temperature of the Tayuan Well in Beijing during multiple earthquakes, scholars Yang et al. [2] believe that in the co-seismic response, the convection and mixing of the water body in the well are the causes of the water temperature drop. Shi et al. [3] studied the co-seismic water level and water temperature changes of the Tangshan Mine Well during 11 earthquakes of magnitude 7.0 or above in 2004, and found that when the well water was vertically oscillated, the disturbed well water triggered a dispersion effect and resulted in changes in the water temperature. During multiple major earthquakes occurred from 2008 to 2013, the water level changes of Jiangsu Well were mostly fluctuations, and the water temperature changes were mostly step changes or trend changes, and the changes of water level and temperature all increased with the increase of the magnitude of the earthquake, and decreased with the increase of the distance from the well to the epicenter [4]. Zhang et al. studied the changes in water level and temperature of Huize Well from 2003 to 2006 and found that there’re obvious relevance, synchronization, and differences between the trend changes of water level and water temperature, and they think such differences are the anomalies before the earthquake [5]. Wang et al. [6] researched the abnormal changes in water level and water temperature and their mechanisms in medium-term, short-term, and short-term stages before the 2014 M4.3 Anhui Shuoshan earthquake. Studies on the fluctuations of water level and water temperature during and before earthquakes play an important role in exploring the formation, occurrence, and development of earthquakes, and such studies can provide important references for pre-earthquake prediction.

Current studies concerning the correlation between water level and water temperature generally use the traditional statistical methods, such as the Pearson's correlation coefficient method, and the linear regression method, etc. However, since the changes in the water level and temperature of monitoring wells are complicated and the time series are unstable, some non-linear correlation methods might be more suitable. At present, the fractal method has been effectively applied to the searching of abnormal soil radon content or electromagnetic anomalies before earthquakes [7-16]. Kar et al. [17] proposed to apply MF-DCCA to study the relationship between soil radon content and temperature anomaly, and they found that temperature anomaly and soil radon content exhibited varying degrees of complexity and multifractal cross-correlation. This paper attempts to use MF-DCCA to explore the non-linear correlation between the water level and water temperature of seismic monitoring wells, namely the multifractal cross-correlation.

In our study, water level and water temperature data of Huize Well from 2004 to 2006 were selected for analysis; at first, the cross-correlation test was carried out to verify that there’s a cross-correlation between water level and water temperature; then, DCCA was used to obtain the cross-correlation coefficient of water level and water temperature under different scales; after that, the MF-DCCA method was adopted to verify the multifractal cross-correlation between the time series of water level and water temperature; at last, this paper used MF-DCCA to figure out the changes in multifractal cross-correlation between water level and water temperature in the sliding window.

3.1 Cross-correlation test

The cross-correlation statistic test can examine whether two groups of time series $\left\{x_{t}, \mathrm{t}=1,2, \cdots, \mathrm{N}\right\}$ and $\left\{y_{t}, \mathrm{t}=1,2, \cdots, \mathrm{N}\right\}$ have long-range cross-correlation [19], wherein the cross-correlation function C_i and the cross-correlation statistic $Q_{c c}(m)$ are defined as follows:

$C_{i}=\frac{\sum_{k=i+1}^{N} \quad x_{k} y_{k-i}}{\sqrt{\sum_{k=1}^{N} \quad x_{k}^{2} \sum_{k=1}^{N} \quad y_{k}^{2}}}$,             (1)

$Q_{c c}(m)=N^{2} \sum_{i=1}^{m} \frac{c_{i}^{2}}{N-i}$             (2)

The statistic $Q_{c c}(m)$ approximately obeys the $\chi^{2}(m)$ distribution with m degrees of freedom. Within the value range of freedom degree m, if statistic  is consistent with the critical value $\chi^{2}(m)$ at the significance level of 5%, then there is no cross-correlation between the two time series; if statistic $Q_{c c}(m)$ is greater than critical value $\chi^{2}(m)$ at the significance level of 5%, then there is a cross-correlation between the two time series.

3.2 DCCA

The DCCA (detrended cross-correlation analysis) method [20] has extended the traditional correlation analysis methods which can only process stationary time series, DCCA can process non-stationary time series, and its specific steps are:

(1) Assume there’re two time series $\left\{x_{t}, \mathrm{t}=1,2, \cdots, \mathrm{N}\right\}$ and $\left\{y_{t}, \mathrm{t}=1,2, \cdots, \mathrm{N}\right\}$, calculate their sums as:

$\left\{ \begin{matrix}   {{R}_{k}}=\sum\limits_{i=1}^{k}{{{x}_{i}},k=1,\cdots ,N}  \\   {{{{R}'}}_{k}}=\sum\limits_{i=1}^{k}{{{y}_{i}},k=1,\cdots ,N}  \\ \end{matrix} \right.$               (3)

(2) Divide the two cumulative time series into N-n overlapping time slots, and each time slot contains n+1 values. The time slot starts at i and ends at i+n.

(3) In each time slot, the two cumulative time series are respectively fitted using the least squares method to obtain trends of local time slots $\widetilde{R_{k, l}}$ and $\widetilde{R_{k, l}^{\prime}}$.

(4) In each time slot, get the residuals of two original cumulative time series and their respective local fitted time series trends, and calculate the covariance of their residuals:

$f_{DCCA}^{2}(n,i)=\frac{1}{n-1}\sum\limits_{k=i}^{i+n}{({{R}_{k}}-\widetilde{{{R}_{k,i}}})}({{{R}'}_{k}}-\widetilde{{{{{R}'}}_{k,i}}})$         (4)

(5) Calculate the covariance of the entire time series on different time scales:

$F_{DCCA}^{2}=\frac{1}{N\text{-}n}\sum\limits_{i=1}^{N-n}{f_{DCCA}^{2}}(n,i)$         (5)

Calculate Formula (5) repeatedly on each time scale, then we can get:

${{F}_{DCCA}}(n)\sim {{n}^{\lambda }}$       (6)

$\lambda$ is the scale parameter, it can reflect the cross-correlation of two time series. When $0<\lambda<0.5$, it means that the two time series have long-range negative correlation; when $0.5<\lambda<1$, it means that the two time series have long-range positive correlation. When $R_{k}=R_{k}^{\prime}$, the detrended covariance $F_{D C C A}^{2}(n)$ becomes the detrended variance $F_{D F A}^{2}(n)$. Since the DCCA scale parameter $\lambda$ cannot strictly quantify the correlation between the two time series, so the detrended variance DFA was added to obtain the detrended correlation coefficient $\rho_{D C C A}$, which is:

$\rho_{D C C A}=\frac{F_{D C C A}^{2}\quad(n)}{F_{X D F A}\quad(n) F_{Y D F A}\quad(n)}$            (7)

where, $F_{X D F A}(n)$ and $F_{Y D F A}(n)$ are the DFA scale parameters of the two time series $\left\{x_{t}, \mathrm{t}=1,2, \cdots, \mathrm{N}\right\}$ and $\left\{y_{t}, \mathrm{t}=1,2, \cdots, \mathrm{N}\right\}$. The DCCA method is different from the traditional Pearson correlation coefficient method, Pearson correlation coefficient is applicable to stationary time series, while the DCCA has extended to non-stationary time series. The traditional Pearson correlation coefficient has a unique value, while DCCA obtains different correlation coefficient values for different time scales.

3.3 MF-DCCA

Zhou [21] combined multifractal method with DCCA and proposed the MF-DCCA method, namely multifractal detrended cross-correlation analysis, it is also a research method for non-stationary time series.

The specific steps are as follows:

(1) For two given time series $\left\{x_{t}, \mathrm{t}=1,2, \cdots, \mathrm{N}\right\}$ and $\left\{y_{t}, \mathrm{t}=1,2, \cdots, \mathrm{N}\right\}$, N is the length of the time series, calculate the cumulative deviation sums of their mean values, and generate new series $X(i)$ and $Y(i)$:

$X(i)=\sum_{t=1}^{i}(x(t)-\bar{x}) \quad, Y(i)=\sum_{t=1}^{i}(y(t)-\bar{y})$             (8)

where, $\bar{x}$ and $\bar{y}$ are the mean values of the two time series.

(2) Divide the two time series $X(i)$ and $Y(i)$ into $N_{S}$ sub-intervals. The length of each interval is s, and there’s no overlap between each interval, wherein $N_{s}=\operatorname{int}(N / S)$. Since for most time series, the length N might not be divided exactly by s, that is, there might be residual data that is less than length s. In order to make full use of the data, the first and last data of the two series $X(i)$ and $Y(i)$ are reversed, and the series are re-divided for more than once, then $2 N_{s}$ subintervals are obtained:

(3) In each sub-interval, the time series are fitted by the least squares method to generate two fitted series $x_{\lambda}$ and $y_{\lambda}$, which are then detrended to obtain the covariance function as follows:

when $\lambda=1,2, \cdots, N_{s}$,

$\begin{align}  & {{F}^{2}}(s,\lambda )=\frac{1}{s}\sum\limits_{i=1}^{s}{(X((\lambda -1)s+i)-{{x}_{\lambda }}(i))} \\ & \times (Y((\lambda -1)s+i)-{{y}_{\lambda }}(i)) \\  \end{align}$                     (9)

when $\lambda=N_{s}+1, \cdots, 2 N_{s}$,

$\begin{align}  & {{F}^{2}}(s,\lambda )=\frac{1}{s}\sum\limits_{i=1}^{s}{(x(N-(\lambda -{{N}_{s}})s+i)-{{x}_{\lambda }}(i))} \\ & \times (Y(N-(\lambda -{{N}_{s}})s+i)-{{y}_{\lambda }}(i)) \\  \end{align}$                   (10)

(4) Take the mean value of the covariance in each sub-interval, and calculate the q-order fluctuation function:

${{F}_{q}}(s)={{(\frac{1}{2{{N}_{s}}}\sum\limits_{\lambda =1}^{2{{N}_{s}}}{({{F}^{2}}{{(s,\lambda )}^{\frac{q}{2}}})})}^{\frac{1}{q}}}$                (11)

when $q=0$,

${{F}_{0}}(s)=\exp (\frac{1}{4{{N}_{s}}}\sum\limits_{\lambda }^{2{{N}_{s}}}{(\ln {{F}^{2}}(s,\lambda ))})$               (12)

(5) The power-law relationship between the fluctuation function $F_{q}(s)$ and time interval s under different scales s is:

$F_{q}(s) \sim s^{H(q)}$, namely: $\log F_{q}(s)=H(q) \log (s)+\log C$

where, $H(q)$ is the generalized Hurst exponent.

(6) The multifractal scale index is $\tau(q)$, and its relationship with the generalized Hurst exponent is:

$\tau (q)=qH(q)-1$

If $\tau(q)$ and $q$ have a non-linear relationship, then the two time series have multifractal characteristics, otherwise they are single fractal.

(7) Through Legendre transformation, the multifractal spectrum is obtained as:

$\alpha (q)=H(q)+{H}'(q)$

$f(\alpha )=q(\alpha (q)-H(q))+1$

where, $\alpha$ is the singularity index of the time series, $f(\alpha)$ is the multifractal spectrum, wherein the width of the multifractal spectrum is $\Delta \alpha=\alpha \min _{\max }$. The greater the width of the multifractal spectrum, the stronger the multifractals of the cross-correlation of the two time series, and the more obvious the long-range cross-correlation; otherwise, the weaker the multifractals and the long-range cross-correlation. $\Delta f=f\left(\alpha_{\min }\right)-f\left(\alpha_{\max }\right)$ is the difference of the multifractal spectra. When $\Delta f>0$, the multifractal spectrum shows left skew, it means that it is of higher probability that the two time series are of greater values, and they have a rising trend; when $\Delta f<0$, the multifractal spectrum shows right skew, it means that that it is of higher probability that the two time series are of smaller values, and they have a decline trend.

Since the high-frequency information of well water level and temperature may contain abnormal information of the earthquakes, and the water level anomalies might be caused by non-tectonic factors, therefore, whether other factors would cause the width of the well water level and temperature multifractal spectrum before the Huize earthquake to increase should be considered. At first, by looking up the changes of the air pressure coefficient, whether the anomaly of the multifractal spectrum width before the Huize earthquake was caused by air pressure could be determined. Based on the air pressure data (in days) from 2004 to 2006, a one-variable linear regression model with a time window of 2 months describing the relationship between air pressure and the water level of Huize well was established to calculate the air pressure coefficient in the response time window.

Table 3. Air pressure coefficient of Huize well

The linear regression models that past the significance test in different time windows, namely in those there’s a regression relationship between the air pressure change and the water level change, are shown in the Table 3. In terms of the air pressure coefficient, before the Huize earthquake, there’s no obvious anomaly in the air pressure coefficient, and its value was between 2.3179 and 3.0521, indicating that the high-value anomaly of the multifractal spectrum width before the Huize earthquake had nothing to do with air pressure; in addition, the correlation between the general atmospheric pressure change and the water level change was poor, only the correlation coefficients of 2004/11/01-2004/12/31 and 2006/05/01-2005/06/30 time intervals were above 0.7, therefore, the impact of air pressure on the water level was average.

6.png

Figure 6. Data of theoretical of earth tide (in hours) of Huize well from 2004/01/01 to 2006/12/31

In order to determine the relationship between the change of earth tide and the change of water level, the data of the theoretical value of earth tide (in hours) from 2004/01/01 to 2006/12/31 were generated as shown in the following figure 6, again, 2 months were taken as a time window, within which, a one-variable linear regression model with earth tide change as the independent variable and the water level change as the dependent variable was established. However, within 18 time windows, the linear model failed to pass the significance test, indicating that the impact of earth tide on water level was not significant.

The rainy season in Huize area is from May to August (Figure 7). Observation of the water temperature of Huize well shows that the well water temperature has a long-term downward trend, which is a normal dynamic change of the water temperature, and it isn’t affected by the rainfall. The pressure well is considered to be well sealed, that is, the situation that the rainfall would cause high-value anomaly to the width of the multifractal spectrum before Huize earthquake has been excluded.  

7.png

Figure 7. Data of rainfall (in days) from 2004/01/01 to 2006/12/31

As a result, possible causes of the anomaly in the width of the multifractal spectrum, such as air pressure, earth tide, and rainfall had all be excluded, therefore, this paper believes that the possible reason for such situation might be that, when the mechanical processes generated by the earthquake spread to the aquifer rock mass in the shallow crust, due to the existence of voids (such as pores, cracks, and caves) developed in the rock mass, even a weak mechanical force can cause certain deformation to the voids, or cause some weak parts to generate local cracks, and such deformation and cracks will certainly lead to changes in the pore pressure of the aquifer, thereby affecting the changes in the water flow inside the aquifer system [22-25]. Since the airtight monitoring wells can sensitively respond to small changes in water level and water temperature, and the changes in the interactive relationship of the fluctuations of well water level and temperature are the critical phenomenon of the groundwater level and temperature before earthquakes, the MF-DCCA method can well capture such fluctuation changes. Before and after the M5.3 Huize earthquake, obvious changes in the width of the multifractal spectrum could be observed, and we believe that this might be the precursor of the Huize earthquake. However, this method needs to be further verified by the changes of water level and temperature of more monitoring wells before the earthquakes.