## Abstract

Due to the influences of geological structures, many confined ascending springs occur in North China's karst area. The electrical conductivity (EC) of karst spring flow is a fundamental variable in characterizing karst systems. However, deeply exploring hidden nonlinear dynamic characteristics is challenging. To avoid overreliance on the fitting polynomial order in the detrending process via classic multifractal detrended fluctuation analysis (MFDFA), the intrinsic time-scale decomposition (ITD) method was applied to identify the external data trend term by decomposing the original data into different frequency modes, and the ITD-MFFA method was proposed to reveal the formation mechanism of spring EC complex characteristics in Jinan's spring area. The results showed that the EC sequences of North China's karst springs are characterized by multifractal behavior and anti-persistence and exhibit multiyear complexity with an overall decreasing trend. Different recharge sources, formation conditions, and seasonal precipitation could be the primary factors driving spring EC complexity. Compared to those of the traditional MFDFA method, the ITD-MFFA method achieved improved anti-interference ability and stability. Timely determination of spring EC data dynamic complexity and monitoring future spring water quality trends can provide guiding significance for protecting karst springs.

## HIGHLIGHTS

Compared with the MFDFA method, the proposed new approach (ITD–MFFA) improves anti-interference ability and stability.

Multifractal fluctuation feature in karst springs’ electrical conductivity data is found.

The complexity of springs' electrical conductivity is not only affected by many factors but also significantly higher in the dry season than in the wet season.

## INTRODUCTION

Due to the influences of geological structures, North China's karst area hosts many confined ascending springs, of which there are mainly two types: igneous rock intrusion contact ascending spring and structural erosion skylight overflow spring. Various hydrogeological factors affect spring water dynamics, which are characterized by nonlinear and complex fluctuations. In recent years, the vulnerability and spring water sensitivity of karst systems have increased under the influence of anthropogenic activities. Consequently, many karst springs present hydrogeological problems such as water quality deterioration and flow attenuation. Due to their high anisotropy and heterogeneity, karstic aquifers contain many voids, such as pores, fractures, and dissolved pores, and a small network of flow channels. Water flowing in such a karstic aquifer exhibits coexisting fissure and channel flows, laminar and turbulent flows, linear and nonlinear flows, and continuous flows and isolated water bodies (Luo *et al.* 2005). Because of karstic systems' complexity and diverse water flow modes, it is challenging to quantify the evolution of their water cycle. Spring water is the natural discharge outlet of karstic aquifers and contains a large amount of geological and hydrological information, manifesting as the equilibrium among various elements within these systems (Manga 2001). As a typical nonstationary time series, a spring's electrical conductivity (EC) curve has intense volatility and variability, showing random chaotic characteristics (Zhao 2018). The nonstationary changes and uncertainty in the EC time series are defined as the EC complexity, and spring EC complexity is much higher than that of other spring water dynamic parameters such as water level and spring flow. External factors such as atmospheric precipitation and artificial source inputs and water–rock interactions can disturb a karstic system's spring water EC, causing its change rate to lag behind water level changes. Hence, spring EC's nonstationarity is a comprehensive characteristic of spring water factors, such as the water level, water quality, flow rate, and water temperature, during the karstic water cycle. Unfortunately, nonstationary time series, such as precipitation variations, suffer from significant prediction and risk management issues considering their nonlinearity and uncertainty (Wen *et al.* 2017). Accordingly, conducting an in-depth investigation of the inherent complex fluctuations of spring EC and its mechanism has a significant theoretical value for revealing karstic systems' nonlinear dynamic evolution.

Previous spring EC series studies were based on linear dynamics statistical methods, such as regression analysis (Meng *et al.* 2021). However, these classic methods can only describe small-scale data's overall trend and mutation characteristics; as the data scale expands, the sequence these methods describe tends to be smooth (Liu 2019). Therefore, it is challenging to describe karstic systems' complex hydrodynamic uncertainty characteristics systems using these methods. The application of complexity science in time series provides vital theoretical and methodological support for solving related problems, deepening the study of groundwater complexity. For instance, Liang *et al.* (2016) systematically determined the fractal characteristics of groundwater-level fluctuations and accurately described the response of those fluctuations to groundwater flow system dynamics. These results also provided new insights for research on spring EC complexity.

To obtain a complete description of complex time series, Kantelhardt *et al.* (2002) proposed the multifractal detrended fluctuation analysis (MFDFA) method based on the detrended fluctuation analysis (DFA). The trend term's interference in time series was eliminated through the detrending process, and various wave functions were considered to describe the fractal structure of time series with different dimensions (Hou *et al.* 2010). From a dynamic perspective, the MFDFA method filters the original sequence's unknown trend components, while the remaining deviation sequence retains its fluctuation components of the original sequence, which has been applied to many disciplines (Zhang *et al.* 2019; Balkissoon *et al.* 2020; Chanda *et al.* 2020; Sarker & Mali 2021). However, the MFDFA method faces certain shortcomings in detrending. For example, in the MFDFA method's detrended fitting process, the high-order trend terms fitted via the least-squares method produce significant errors. Furthermore, the original sequence's cumulative difference processing frequently causes valuable information to be lost. Effectively removing the external trend items in the original data sequence is the main difficulty and challenge. Therefore, the data decomposition method is considered to decompose the original data into multiple modal components and perform data reorganization to solve the shortcomings of the traditional MFDFA method in the detrending process. The commonly used data decomposition methods include wavelet transform (WT), empirical mode decomposition (EMD), and variational mode decomposition (VMD). The WT method is based on the original signal for the composite resolution analysis, but it cannot separate similar frequencies and lacks adaptability (Wang *et al.* 2018); although the EMD method can adaptively decompose the nonstationary time series into multiple intrinsic mode components, this method still has shortcomings such as endpoint effect and mode aliasing, which significantly affect the performance of the reconstructed data. The VMD method transforms the modal decomposition problem into a constrained variational problem, which solves the shortcomings of the EMD method (Lin *et al.* 2022). However, the VMD results depend on multiple input parameters (such as the number of decomposition modes), and the arbitrary determination of input parameters will introduce errors into the VMD. Alternatively, intrinsic time-scale decomposition (ITD) (Frei & Osorio 2007) is an adaptive time–frequency analytical method that is superior to the aforementioned methods in suppressing endpoint effects and computational efficiency. This method does not require envelope fitting iteration, and only linear interpolation is performed, so the operation is faster, and the modal mixing is improved compared to EMD and its extended algorithm (Ma *et al.* 2021). This method has been successfully applied to data preprocessing in signal recognition and mechanical fault diagnosis, but its application in the complexity measurement of karst water systems is rarely involved.

To further improve the accuracy and reliability of complex feature recognition and avoid the instability of the MFDFA method in the process of detrending, the ITD method is introduced to identify the external trend term of the data by decomposing the original data into different frequencies. The specific objectives of this paper include (1) proposing the ITD–MFFA method for accurately identifying and quantifying spring EC complexity; (2) quantifying the multifractal characterization of spring EC data in Jinan, China; (3) analyzing the spatiotemporal variability in the complexity of spring EC data and exploring its possible environmental causes; and (4) presenting the major conclusions drawn from the study.

## STUDY AREA

*et al.*2021). The geological structure is monoclinal, with Paleozoic Cambrian and Ordovician carbonate strata as the main body. In the north, the Quaternary strata in front of the mountains conceal the Ordovician strata that connect with the Yanshanian igneous rocks (Zhu

*et al.*2017). After atmospheric precipitation recharge in the southern limestone distribution area, karst water migrated from south to north along the strata and was blocked and enriched by magmatic rocks. The karst water was exposed in the form of springs in the suitable terrain and favorable structural parts (Xing

*et al.*2017). There are 118 springs in an area of 2.6 km

^{2}, which are divided into the Baotu, Zhenzhu, Wulongtan, and Heihu spring groups (Zhu

*et al.*2020). The principal aquifer comprises Cambrian dolomite and dolomitic limestone. Under natural conditions, spring water is the primary drainage form of the karstic system. From 1972 to 2003, spring water was intermittently cut off due to excessive groundwater exploitation. Since 2003, measures such as artificial replenishment and groundwater source closure have been implemented, and artificial source replenishment during the dry season has been conducted in the limestone of Zhangxia Formation underlying the Yufu River and the Quanlu River with time (Xing

*et al.*2018).

## MATERIALS AND METHODS

### Data collection

Considering the spring water exposure position and the difficulty of sampling, this paper takes the Baotu spring, Heihu spring, Wangfuchizi, and Tanxi spring as representatives of the Baotu, Heihu, Zhenzhu, and Wulongtan spring groups, respectively, to obtain conductivity data. The spring water level data were recorded by a wireless telemetry water level meter at intervals of 1 day with an accuracy of ±1 mm. The spring water EC was measured once a week in the field using an *in situ* Aqua TROLL600 multiparameter detector made in the USA, and the EC precision was ±0.1 μS cm^{−1}. After the measurement, through the VuSit™ mobile phone application, monitoring logs and data files were established and summarized in the mobile device to facilitate later download. Due to instrument failure and other reasons, to ensure the continuity of data, this paper collected and collated water level and conductivity data from 2013 to 2020 as the research object. The daily precipitation data (2013–2020) were obtained from 19 rain stations (Figure 1).

### ITD method

*et al.*2021). An operator

*β*is defined as the baseline extraction operator, when

*β*acts on the original sequence, the original sequence is decomposed into low-frequency and high-frequency signals:where

*L*represents the low-frequency component of the original signal, and

_{t}*H*is the proper rotation (PR) component.

_{t}- 1.
- 2.
*β*_{t}is taken as the original signal, and the aforementioned steps are repeated for decomposition until the baseline component is monotonic, where the decomposition is stopped. - 3.

### MFDFA method

The MFDFA method can be used to determine whether multifractal characteristics exist in time series (Zhang *et al.* 2019). For a time series *x*(*t*)*, t* = 1, 2, … , *N*, where *N* is the time-series length. The calculation process is as follows (Kerkhoven & Gan 2011):

*Step 2*: Starting from the sequence's head and tail {

*y*(

*t*)}, two isometric divisions of length

*S*are conducted to obtain 2

*Ns*disjoint subsequences, where

*N*s = int(

*N*/S). If the number of nonoverlapping subsequences is defined as

*v*(

*v*= 1, 2, … , 2

*N*

_{S}), the elements of subsequence

*v*are recorded as {

*y*(I

*+*

_{v}*t*),

*t*= 1, 2,}. Trend function

*y*(

_{v}*t*) for each subsequence can be obtained by polynomial fitting (least-squares regression):where

*b*

_{0},

*b*

_{1},

*b*

_{2},

*b*

_{3}, … ,

*b*

_{m}are the polynomial's coefficients and

*m*is fitting polynomial's order.

If *h*(*q*) is a linear function of *q*, the time series has monofractal characteristics. If *h*(*q*) is a convex function of *q*, the image of log_{2}*F*(*s*)–log_{2}(*S*) is a group of straight lines with unequal slopes, and the time series has multifractal characteristics (Feng 2018).

### Multifractal strength analysis

*h*(

*q*) and

*τ*(

*q*) can be obtained by using Equation (11) (Koscielny-Bunde

*et al.*2006), which is as follows:

If *τ*(*q*) and *q* are linearly related, the time-series data are monofractal. However, if *τ*(*q*) and *q* are nonlinearly related, the time-series data are multifractal, and the stronger their nonlinear relationship is, the stronger are the multifractal characteristics.

The meaning of *α* is as follows: if *α* = 1, the time-series data distribution is uniform; if *α* > 1, the singularity degree is small; if *α* < 1, the singularity degree is large.

*α*

_{max}and

*α*

_{min}are the maximum and minimum of

*α*, respectively, and Δ

*α*denotes the uneven distribution of the measured probability on the overall fractal structure. Specifically, the larger the value of Δα is, the more disordered the time series distribution and the stronger its multifractal features are.

### Modified ITD–MFFA method

Based on the traditional MFDFA method, polynomial fitting is applied to detrend the time series. However, during polynomial fitting, valuable information of the original sequence is inevitably lost. Therefore, the ITD–MFFA method proposed in this study is employed using the following specified steps:

First, the ITD method is adopted to identify and decompose the original sequence and remove the low-frequency signal to form a new time series.

Second, because the external low-frequency trend term is deleted in the first step, the detrending process does not need to be repeated when calculating the multifractal step. Therefore, Equation (6) is omitted, *y _{v}*(

*i*) in Equations (7) and (8) is eliminated, and the original sequence's MFDFA process is simplified to the MFFA process.

### Selecting the parameters for calculating the complexity of spring EC fluctuations

In the time series, the trend of spring EC is judged by log_{2}*F*(s)–log_{2}(*S*)curve. The slope of the curve is *h*(*q*), where *h*(2) has an important analysis significance: when 0.5 < *h*(2) < 1, the sequence shows an overall growth trend, and when the *h*(2) is greater, the sequence growth trend is stronger; when 0 < *h*(2) < 0.5, the sequence showed a downward trend, when the *h*(2) is smaller, the sequence downward trend is more obvious; in particular, when *h*(2) = 0.5, the sequences are random. If the *h*(*q*) and multifractal mass index *τ*(*q*) show nonlinear variations with *q*, the EC time series exhibits multifractal characteristics. The complexity of conductivity change is measured by the width of the multifractal spectrum (Δ*α*).

Currently, there is no clear definition of the length of the nonoverlapping subinterval *S* and the value of *q*. The choice of *S* and *q* in this paper still draws on previous research results. It is generally considered that the minimum value of *S* should be between 10 and *N*/4 of the sequence length (Zhao 2018). Therefore, to obtain a highly stable *F*(*s*), considering the measured sequence length in this study, the *S* length value is limited to [10, 100].

For the general time-series study, the maximum value of |*q*| is 5 (Ihlen 2012). Therefore, this study determined the multifractal order *q* value between −5 and 5, crossing the 0-point at the center to comprehensively investigate the impacts of large and small sequence fluctuation deviations on the wave function and sequence's multifractal characteristics.

## RESULTS AND ANALYSIS

### Dynamic nonlinear characteristics of long spring EC time series

^{−1}).

*et al.*2019). The parameter test method knows the distribution of the data to be tested in advance and tests whether the distribution parameters fall within the corresponding threshold range. The nonparametric test method uses the sample data distribution to test the distribution of the overall data when the data distribution is unknown (Gong

*et al.*2019). The bispectrum analysis selected in this paper is a kind of nonparametric test method. The detection process does not require prior knowledge of the process and only needs to collect process data. Compared with the parametric test method, the former has the advantages of simple operation, fast calculation speed, and wide application range. According to the previous studies, the bispectrum of a time-series subject to a normal distribution is relatively uniform and horizontal, and there are no notable fluctuations, and a random normal bispectrum fluctuates around 0. The bispectral analysis of random nonlinear data yields uneven side views, large fluctuations, and clear hierarchal features, whereas the top view exhibits obvious symmetry (Yu 2013). From Figure 4, although the results of the bispectral analysis of each spring group's spring water EC time series are not the same, they are all nonflat spectra, and their top views reveal obvious symmetry. Therefore, the spring water EC data of the four spring groups are nonlinear time series with typical complex characteristics, which can be considered to analyze the system's complexity further.

### Analysis of the dynamic complexity of the spring EC

_{2}

*F*(

*s*)–log

_{2}(

*S*) curves of the mutually independent subinterval

*S*and

*q*-order wave function

*F*(s) obtained from Equations (9) and (10), respectively, comprise a group of straight lines with different slopes. With the increasing

*q*, the fitting-line spacing decreases, indicating that the generalized Hurst index exhibits a convergent trend with an increasing order

*q*. Moreover, the spring EC time-series' multiscale relationship features and multifractal characteristics can be observed (Figure 5(a)), and nonlinear curves of

*h*(

*q*) and

*τ*(

*q*) with the

*q*value are obtained from Equations (11)–(13). For

*q*= 2,

*h*(

*q*) < 0.5 (Figure 5(b) and 5(c)) means the sequence achieves antipersistence and a long memory in the long-term fluctuation process. Changes in historical data, current data, or future changes are not isolated. The past data fluctuation affects the current data, and this memory might persist for long. Figure 5(d) shows that the EC time-series' multifractal spectrum exhibits a significant single-peak bell shape, and the spectrum's peak is to the right of the center, revealing the decreasing trend of the spring EC time series at this stage. The spectrum Δ

*α*width is 1.6095, and the nonuniformity of the data series' fractal structure is high, indicating that the data fluctuate notably with a high degree of complexity and notable multifractality (Zhao 2018).

When affected by atmospheric precipitation, Jinan spring water has the characteristics of short-term concentrated recharge in the wet season and long-term consumption in the dry season (Xing *et al.* 2018). Previous studies have shown that spring levels are correlated to varying degrees with precipitation in the previous three hydrological years (Meng *et al.* 2021). Theoretically, the greater the amount of data used, the more realistic the complexity analysis results of the karst spring system. The short time period selected during the analysis process will result in insufficient information. To fully explore the evolution of spring EC multifractal characteristics over the past 8 years (2013–2020), every four consecutive hydrological years were selected as a basic sequence unit, and the semiannual sliding recursive method was used to divide the original time series into the following seven periods: 2013.6–2017.6, 2014.1–2018.1, 2014.6–2018.6, 2015.1–2019.1, 2015.6–2019.6, 2016.1–2020.1, and 2016.6–2020.6. Because the information contained in the *h*(*q*)–*q* relationship curve and multifractal spectrum is comprehensive, it can be directly applied to determine whether the spring EC time series exhibits complex fractal characteristics. In this study, the *h*(*q*)–*q* relationship curve and multifractal spectrum are directly generated at various stages (Figure 5(e)–5(f)). During these periods, the variation in *h*(*q*) with *q* reveals obvious nonlinearity, and the multifractal spectrum demonstrates symmetric, single-peak bell curve characteristics, indicating that the spring EC time series during each period display complex multifractal characteristics. Furthermore, the spring EC Hurst exponent during each period is <0.5 for *q* = 2, indicating that the series during each period exhibit antipersistence. Therefore, the future increment in the system is negatively correlated with the past increment.

In the multifractal analysis, the singular index *α* reflects the fluctuation phenomenon's local irregularity in the sequence. The larger the value is, the lower is the singularity and the smoother is the data sequence in this part. Δ*α* is also referred to as the multifractal characteristic parameter or fingerprint (Zhao 2018).

Spring . | Max. (μS·cm^{−1})
. | Min. (μS·cm^{−1})
. | Mean (μS·cm^{−1})
. |
---|---|---|---|

Heihu spring | 954 | 831 | 895 |

Baotu spring | 858 | 722 | 777 |

Tanxi spring | 776 | 671 | 707 |

Wangfuchizi | 718 | 622 | 685 |

Spring . | Max. (μS·cm^{−1})
. | Min. (μS·cm^{−1})
. | Mean (μS·cm^{−1})
. |
---|---|---|---|

Heihu spring | 954 | 831 | 895 |

Baotu spring | 858 | 722 | 777 |

Tanxi spring | 776 | 671 | 707 |

Wangfuchizi | 718 | 622 | 685 |

To further examine the fluctuation in the complexity of the EC time series with time, the calculation results of the multifractal spectrum of the EC time series during different periods were used to calculate the statistics of the complexity index Δ*α*, and the Δ*α* curve was generated (Table 2). The spring EC complexity indicates a decreasing trend with time; however, the complexity exhibits some fluctuations during the different periods. In the periodic change in the spring water EC of the Baotu spring, the complexity of the spring water EC reveals an obvious decreasing trend from June 2017 to January 2019 and from June 2019 to June 2020 rather than from June 2013 to January 2015 and from June 2015 to June 2016. The spring water EC complexity exhibits an increasing trend from January 2019 to June 2019 compared with January 2015 to June 2015. The fluctuation state of the complexity of the EC time series coincides with the spring's antipersistence, whereas the EC's Hurst index during the different periods is <0.5.

Period . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . |
---|---|---|---|---|---|---|---|

Δα | 1.42 | 1.35 | 1.34 | 1.13 | 1.29 | 1.03 | 1.06 |

Period . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . |
---|---|---|---|---|---|---|---|

Δα | 1.42 | 1.35 | 1.34 | 1.13 | 1.29 | 1.03 | 1.06 |

### Spatial analysis of the dynamic complexity of the spring EC

*h*(

*q*) decreases with the increasing

*q*; all four spring groups' spring EC time series exhibit multiscale and multifractal characteristics. The

*h*(2) value reflects the long-term correlation strength. All corresponding

*h*(2) values of the spring groups are <0.5, indicating that each spring group's EC time series displays a negative long-term correlation, i.e., antipersistence. The

*h*(2) values decrease in the order of Heihu spring > Tanxi spring > Wangfuchizi spring > Baotu spring; in addition, the antipersistence is the highest for Baotu spring's EC. Figure 6(b) and Table 3 show the multifractal spectrum of the EC time series and the singular spectrum width Δ

*α*for each spring group. The variation range is 1.56–1.81, indicating that each spring group's multifractal degree varies. The Wangfuchizi spring's Δ

*α*value is the largest, indicating that the corresponding EC time series attain the highest multifractality degree, whereas Tanxi, Baotu, and Heihu springs' values considerably fluctuate.

Spring . | Long-term correlation exponent h(2)
. | Complexity Δα
. |
---|---|---|

Baotu spring | 0.316 | 1.606 |

Heihu spring | 0.331 | 1.559 |

Tanxi spring | 0.374 | 1.722 |

Wangfuchizi | 0.3839 | 1.814 |

Spring . | Long-term correlation exponent h(2)
. | Complexity Δα
. |
---|---|---|

Baotu spring | 0.316 | 1.606 |

Heihu spring | 0.331 | 1.559 |

Tanxi spring | 0.374 | 1.722 |

Wangfuchizi | 0.3839 | 1.814 |

## DISCUSSION

### Interaction mechanism among the potential factors of the spring EC complexity

According to the aforementioned complexity calculation results, spatially, the complexity is low in the south (Baotu and Heihu springs) and high in the north (Tanxi and Wangfuchizi springs). Different zones and influencing factors in this complex and dynamic groundwater system impose varying effects on the karstic system. Atmospheric precipitation, the formation conditions of spring water, and the spring source might affect the spring water's dynamic complexity. These factors interact and restrict each other, forming the basic pattern of spring EC complexity. Therefore, based on the preliminarily determined spring EC complexity difference, the study examined the effects of atmospheric precipitation, the spring formation conditions, and the spring source on the variations in spring EC complexity.

*et al.*2015; Rajaee

*et al.*2019; Pandey

*et al.*2020; Afan

*et al.*2021; Band

*et al.*2021), the water level can be accurately predicted. Although the conductivity data of the spring water in the diagram show the seasonal periodic variation characteristics of increasing in the wet season and decreasing in the dry season similar to the fluctuation of the water level, it only fluctuates up and down within a certain range, and there is no strict correlation with the water level. Atmospheric precipitation is an important factor in the complex fluctuation of the conductivity of spring water, but it is not the only factor. The fluctuation of the conductivity of the spring water is more complex than the fluctuation of the water level. To fully explore the influence of atmospheric precipitation in different seasons on the complexity of spring water conductivity, the complexity of spring water conductivity in wet and dry seasons is calculated (Figure 7(b)). The results showed that the complexity of spring water conductivity in the dry season is significantly higher than that in the wet season, suggesting that the internal effect of the spring water dynamic system is more complicated in the dry season. Therefore, it is necessary to further explore other influencing factors of the complex fluctuation of spring water conductivity in combination with hydrogeological factors such as spring water supply source and formation conditions.

The difference of recharge source and water circulation depth caused the fluctuation of spring conductivity complexity. The Jinan Spring Group is jointly supplied by the Ordovician-Fengshan Formation aquifer and the Zhangxia Formation aquifer (Zhu *et al.* 2020). There are differences in the hydrochemical types of groundwater in different strata. For example, the hydrochemical types of groundwater in the Cambrian Zhangxia Formation are mainly HCO_{3}·SO_{4}-Ca and HCO_{3}-Ca, and the average TDS is 498mg/L. The hydrochemical types of groundwater in the Ordovician–Cambrian Fengshan Formation are mainly HCO_{3}·SO_{4}-Ca and HCO_{3}-Ca·Mg, and the average TDS is 512 mg/L (Wu 2021). The karstic water-bearing medium of the Jinan spring primarily comprises dissolution fractures (Li *et al.* 2021). There are also obvious differences in the process of the karst water cycle. Only minor flow channels and karst caves are developed in the Ordovician strata at a burial depth of 100 m (Zhu *et al.* 2020), and the degree of the karst development gradually decreases with the increasing depth. According to the previous research results, the shallow karst development is high with the buried depth of 100 m, and the retention time of the karst water is less than 6 years. The retention time of the karst water in Ordovician–Cambrian Fengshan Formation with the buried depth of more than 100 m is 6–55 years, and the retention time of groundwater in Zhangxia Formation aquifer is 48–82 years (Meng 2021). The shallow karst water cycle is rapid, and the update of water quantity and water quality is fast, and so the concentration of Cl^{−}, , and in the spring water is small. With the increase of runoff time, water–rock interaction and dissolution are intensified (Cao 2020), and the difference of the karst water quality is more obvious. According to the EC measurement results, the spring water EC values vary between the Cambrian Fengshan Formation–Ordovician karst aquifer and the Zhangxia karst aquifer (Hou 2019). The karst water from different strata converges at the spring outlet, causing conductivity fluctuations.

Seasonal variation of the spring runoff channel also has some influence on conductivity fluctuation. According to the seasonal periodic variation characteristics of EC and the groundwater migration path profile (Figure 7(c)), the spring water level rises during the wet season. Spring water is mostly fed from the Cambrian Fengshan Formation–Ordovician aquifer. Large karst fissures at shallow depths are quickly recharged via atmospheric precipitation, thereby collecting and conducting water. The primary spring water channel is a large fissure in the aquifer at this stage. The spring water level is low during the dry season, and the spring water originates primarily from the Cambrian Zhangxia Formation aquifer. Due to the weakening of the vertical karst development, maintaining the spring water level is dependent on the slow release of water from fissures and water-bearing pores in the karst aquifer. Furthermore, to prevent spring water from being cut off, artificial sources in the Zhangxia Formation limestone underlying Yufu and Quanlu Rivers have been supplemented during the dry season with time (Xing *et al.* 2018). However, mixing with this supplementary water intensified the spring EC fluctuation complexity.

The difference of geological origin has a direct influence on the complexity of spring conductivity. In terms of its origin, the Jinan spring's geological structure is a north-dipping monocline dominated by the Paleozoic strata. Spring water flows from the southern mountainous area toward the north and is blocked by the Yanshanian intrusive complex, allowing spring water to accumulate. Due to the heterogeneous spatial distribution of these intrusive rocks (Figure 7(d)), each of the four spring groups' geological sources differs. The intrusive rocks in some regions of the Baotu and Heihu springs are entirely denuded, whereas carbonate rocks are in direct contact with the upper gravel layer; thus, water from the underground aquifer can ascend to the surface to form springs. The diorite overlying the limestone in the Wulongtan (Tanxi spring) and Zhenzhu (Wangfuchizi spring) spring groups serves as a poor aquitard for this spring water. Furthermore, the fault or fracture zone passing through the diorite functions as a channel, allowing water to surge upward from the deep limestone aquifer. According to their different geological origins, the Baotu and Heihu spring groups comprise ascending springs emerging from igneous rocks, whereas the Wulongtan (Tanxi spring) and Zhenzhu (Wangfuchizi spring) spring groups contain overflowing springs due to igneous rock erosion. According to the calculation results of the spring EC complexity, the complexity of the Baotu and Heihu springs is low, whereas that of the Tanxi and Wangfuchizi springs is high. Thus, differences in the geological origin directly impact spring EC complexity.

### Analysis of the performance of the ITD–MFFA algorithm

#### Stability of the ITD–MFFA algorithm

*m*is set to 1, 2, 3, and 4, the spring EC complexity values computed using the MFDFA algorithm are 1.845, 1.831, 1.802, and 1.821, respectively. The results of the MFDFA algorithm depend on the selection of trend term parameters, which causes the instability of the algorithm. The ITD–MFFA method decomposes the data trend of the original data at first to effectively avoid the influence of different parameters on the results. Therefore, the complexity analysis results based on the ITD–MFFA method are more stable than those based on MFDFA.

#### Reliability of the ITD–MFFA algorithm

To evaluate the ITD–MFFA algorithm's reliability, for verification, white noise and colored noise with an intensity of 1 dB and amplitude of 5% of the sequence standard deviation are added to the original sequence without adjusting the original sequence's attributes (Cheng 2018). The ITD–MFFA and MFDFA methods' calculation results before and after the addition of noise to the original data are compared, and the robustness to noise (i.e., the anti-interference ability) is analyzed. The complexity calculation results obtained using the two methods before and after noise addition are presented in Table 4 and plotted in Figure 8(b), showing that under noise interference, both methods' calculation results exhibit deviations; however, the deviation in the calculation results obtained using the MFDFA method is greater. For the ITD–MFFA method, the correlation coefficients of the sequence complexity before and after the addition of noise are 0.9991 and 0.9993, respectively, and for the MFDFA method, they are smaller at 0.9950 and 0.9355, respectively. Consequently, the ITD–MFFA method displays a stronger anti-interference ability, and the calculation results are more stable. However, due to the detrending process's imprecision, the complexity of the spring EC time series is frequently overestimated using the traditional MFDFA method. Furthermore, the ITD–MFFA method improves the stability and reliability of spring EC measurements.

Spring . | Complexity results with ITD–MFFA . | Complexity results with MFDFA (m = 1). | ||||
---|---|---|---|---|---|---|

Original sequence . | Add white noise . | Add colored noise . | Original sequence . | Add white noise . | Add colored noise . | |

Baotu spring | 1.606 | 1.6063 | 1.6061 | 1.8449 | 1.8453 | 1.8103 |

Heihu spring | 1.559 | 1.5671 | 1.5534 | 1.6531 | 1.6326 | 1.6229 |

Tanxi spring | 1.722 | 1.7231 | 1.7275 | 1.8662 | 1.8437 | 1.8891 |

Wangfuchizi | 1.814 | 1.8126 | 1.814 | 1.9226 | 1.9009 | 1.9336 |

Correlation coefficient R with the original sequence complexity | — | 0.9991 | 0.9993 | – | 0.9950 | 0.9355 |

Spring . | Complexity results with ITD–MFFA . | Complexity results with MFDFA (m = 1). | ||||
---|---|---|---|---|---|---|

Original sequence . | Add white noise . | Add colored noise . | Original sequence . | Add white noise . | Add colored noise . | |

Baotu spring | 1.606 | 1.6063 | 1.6061 | 1.8449 | 1.8453 | 1.8103 |

Heihu spring | 1.559 | 1.5671 | 1.5534 | 1.6531 | 1.6326 | 1.6229 |

Tanxi spring | 1.722 | 1.7231 | 1.7275 | 1.8662 | 1.8437 | 1.8891 |

Wangfuchizi | 1.814 | 1.8126 | 1.814 | 1.9226 | 1.9009 | 1.9336 |

Correlation coefficient R with the original sequence complexity | — | 0.9991 | 0.9993 | – | 0.9950 | 0.9355 |

### Guide to the spring EC complexity

Different geological structures and factors influence the spring water EC fluctuations to varying degrees in a complex karst groundwater system. Studies have reported that karst hydrogeological conditions significantly affect the groundwater EC, comprehensively reflecting regional hydrodynamic conditions and lithologic characteristics (Zhu *et al.* 2019). Due to the spatial heterogeneity of geological structures, the development degrees of groundwater runoff channels and karst landscapes considerably vary over a small area. Furthermore, the spring EC variation is complex because of the uncertainty in the recharge mode and source. Differences in the water–rock contact time driven by seasonal precipitation changed the dynamic structure of the underground aquifer system, and the recharge source for spring water varies, while the dominance of large fissure flow channels and pore–fissure flow alternates. Moreover, the recharge source circulation at different depths and the mixing of spring water with artificially supplemented water affect the spring EC complexity.

Studies have used long-term EC monitoring data to mine information on basin-scale groundwater dynamic conditions and water–rock characteristics and to accurately calculate the proportion of water being supplied to karst springs (Guo *et al.* 2018; Zhu *et al.* 2019; Zhu *et al.* 2020). Jiang *et al.* (2008) and others conducted studies on runoff generation and concentration in surface rock zones combined with EC data. The EC can reflect the precipitation recharge level and the influence of recharge precipitation in the surface karst zone. Furthermore, information on the dynamic changes in the EC is critical for studying karst formation mechanisms and dynamic processes.

Hydrochemically, EC is a vital indicator of spring water quality. The magnitude of the EC reflects the concentration of charged ions in water. The lower the EC value is, the lower is the degree of water pollution (Yu *et al.* 2020). In the Jinan spring, the measured time-series analysis of the spring water EC complexity reveals that the spring water EC dynamics attain a long-range correlation and antipersistence. Therefore, the spring water EC should follow the opposite trend in the future, indicating that the spring water quality will develop along the opposite direction, either good or bad. The higher the spring EC complexity is, the more challenging it is to predict the future water quality, the greater is the number of potential influencing factors, the higher is the uncertainty, and the greater is the water quality risk. Therefore, measuring the spring EC complexity could guide water quality monitoring and management, steer the prediction of the spring water quality, and determine the key areas or directions for the future work.

## CONCLUSIONS

With continuous changes in the groundwater of karstic systems, the complex variations in karst springs' EC represent crucial information. Accurately identifying the complexity of the spring conductivity fluctuation process is of great significance for exploring the dynamic characteristics of the karst water system. Due to the detrended process of the traditional MFDFA method incompleteness, the evaluation results often overestimate the spring EC complexity. To improve the accuracy of the calculation results, this paper proposes the ITD–MFFA approach based on the original MFDFA method and analyzes the temporal and spatial distribution characteristics and potential influencing factors of the EC complexity of North China's karst springs. The ITD–MFFA method improves the accuracy and reliability of analyzing the spring EC complexity. The EC fluctuation in the Jinan karst spring exhibits multifractal characteristics and antisustainability, and the complexity variation reveals a nonsignificant decreasing trend with time. Seasonal differences in spring water sources, the formation conditions of springs, and atmospheric precipitation are most probably the primary driving factors of spring water EC complexity. In a highly variable karst environment, determining the dynamic complexity and future change trends of the spring EC could become a reliable means to understand the formation mechanism and evolution of karst springs. These results can provide important reference information for regional groundwater environmental protection and guiding the restoration of karst springs.

Such a model improves the detrending process of the traditional MFDFA method, avoids the excessive dependence of the calculation results on the polynomial fitting order, and improves the accuracy and stability of the spring dynamic complexity analysis results. However, the selection of the sequence segment length subintervals *S* and fitting order *q* still borrows from previous research results. In addition, due to the limited data collection, some factors that may have a deeper impact on the complexity of spring conductivity have not been determined, such as temperature, which can be examined in the future studies.

## ACKNOWLEDGEMENTS

Thanks to Geological Environmental Research Institute, the Shandong Division of China Metallurgical Geology Bureau and Shandong Provincial Geo-mineral Engineering Exploration Institute for their technical and data support for this study.

## FUNDING

This research was funded by the National Natural Science Foundation of China (41772257 and 42272288) and an innovation team project of colleges and institutions in the Shandong province (2018GXRC012).

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

*Study on the Difference Characteristics of Hydrochemical Dynamics of the Four Springs in Jinan*

*Master Thesis*

*Study on Driving Effect of Regional Water Resources System Complexity Characteristics to Drought Risk and Optimal Allocation*

*Master Thesis*

*Carbon Futures Price Forecasting Based on Improved Support Vector Machine Under Multifractal Principle*

*Master Thesis*

*Study on the Sources of Spring Water Supply in Jinan*

*Master Thesis*

*Research and Application of Time-Varying Hurst Parameter Estimation Methods for Time Series of Self-Similarity*

*Master Thesis*

*Study on the Characteristics of Groundwater Flow System in the Direct Recharge Area of Baotu Spring Area in Jinan*

*Master Thesis*

*Study on the Formation of Chemical Composition of Jinan Karst Spring Water*

*Master Thesis*

*Fractal Characteristics Research of Regional Water and Soil Resources System of Jiansanjiang Branch Bureau*

*Master Thesis*

*Multifractality and Nonstationarity Analysis Based on High Frequency Water Level Data*

*PhD Thesis*