An improved multifractal detrended fluctuation analysis method for estimating the dynamic complexity of electrical conductivity of karst springs

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.

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⁻¹. 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).

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):

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₂F(s)–log₂(S) is a group of straight lines with unequal slopes, and the time series has multifractal characteristics (Feng 2018).

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.

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.

In the time series, the trend of spring EC is judged by log₂F(s)–log₂(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.

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).

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.

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.

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₃·SO₄-Ca and HCO₃-Ca, and the average TDS is 498mg/L. The hydrochemical types of groundwater in the Ordovician–Cambrian Fengshan Formation are mainly HCO₃·SO₄-Ca and HCO₃-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.

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.

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.

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.

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.

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).

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

The authors declare there is no conflict.