Fractal Scaling Analysis of Groundwater Dynamics in Confined Aquifers

Groundwater closely interacts with surface water and even climate systems in most hydro-climatic settings. Fractal scaling analysis of groundwater dynamics is of significance in modeling hydrological 10 processes by considering potential temporal long-range dependence and scaling crossovers in the groundwater level fluctuations. In this study, it is demonstrated that the groundwater level fluctuations of confined aquifer wells with long observations exhibit site-specific fractal scaling behavior. Detrended fluctuation analysis (DFA) was utilized to quantify the monofractality; and Multifractal detrended fluctuation analysis (MF-DFA) and Multiscale Multifractal Analysis (MMA) were employed to examine the multifractal 15 behavior. The DFA results indicated that fractals exist in groundwater level time series, and it was shown that the estimated Hurst exponent is closely dependent on the length and specific time interval of the time series. The MF-DFA and MMA analyses showed that different levels of multifractality exist, which may be partially due to a broad probability density distribution with infinite moments. Furthermore, it is demonstrated that the underlying distribution of groundwater level fluctuations exhibits either non-Gaussian 20 characteristics which may be fitted by the Lévy stable distribution or Gaussian characteristics depending on the site characteristics. However, fractional Brownian motion (fBm), which has been identified as an appropriate model to characterize groundwater level fluctuation is Gaussian with finite moments. Therefore, fBm may be inadequate for the description of physical processes with infinite moments, such as the groundwater level fluctuations in this study. It is concluded that there is a need for generalized governing 25 equations of groundwater flow processes, which can model both the long-memory behavior as well as the Brownian finite-memory behavior.


Introduction
Groundwater in both confined and unconfined aquifers is usually a complex and dynamic system that highly interacts with surface water and even climate systems in most hydroclimatic settlings due to its discharge to rivers and streams and its recharge affected by various related physical processes, such as precipitation, evapotranspiration, and infiltration (Green et al., 2011;Joelson et al., 2016;Li and Zhang, 2007;Rakhshandehroo and Amiri, 2012;Taylor et al., 2013).These processes, which take place over various spatiotemporal scales, add further complexity to groundwater systems.
Groundwater level fluctuations dynamically reflect the responses of an aquifer to its diverse inputs and outputs.Consequently, groundwater level fluctuations are often nonstationary, rendering variabilities over different spatial and temporal scales and resulting in no dependence on single representative spatial and temporal scales.Therefore, groundwater level fluctuations are often characterized as scale-free processes and modeled as fractional Brownian motion (Hardstone et al., 2012;Yu et al., 2016).Although not necessarily totally random, groundwater level fluctuations may demonstrate long-range dependence through time, implying Published by Copernicus Publications on behalf of the European Geosciences Union.
T. Tu et al.: Fractal scaling analysis of groundwater dynamics in confined aquifers a power-law relationship over a variety of timescales, which can be represented by fractals (Yu et al., 2016).
Fractal analysis of both persistent and anti-persistent behavior has been extensively utilized to investigate possible relationships in variability among various scales (Blöschl and Sivapalan, 1995).Temporal fractal scaling analysis of groundwater dynamics can be essential to a better understanding of the modeling of hydrological processes by considering temporal correlations and scaling cascading issues, since groundwater closely links to surface water in hydrological modeling and hydrological models are built upon certain temporal and spatial scales (Blöschl and Sivapalan, 1995;Yu et al., 2016).Hence, fractal scaling analysis of groundwater level fluctuations can guide more representative modeling in hydrological models and in coupled land-atmosphere models.In fact, groundwater dynamics were found to provide a positive feedback to the memory of land surface hydrological processes in climate systems, and enhanced knowledge of fractal behavior in subsurface hydrological processes can help improve weather forecast and climate prediction on different temporal scales (Lo and Famiglietti, 2010).Furthermore, fractal scaling analysis of groundwater level fluctuations may help investigate extreme events and anthropogenic forcing in Earth systems (Yu et al., 2016).
Detrended fluctuation analysis (DFA), originally used to analyze the long-range power-law correlations (i.e., persistent fractal scaling behavior) of time series, is considered a powerful method to quantify the scaling parameter or the Hurst exponent for its capacity in detecting nonstationarities and distinguishing seasonal oscillations from intrinsic fluctuations compared with conventional methods, such as R/S analysis or the variation method (Dubuc et al., 1989;Hardstone et al., 2012;Shang and Kamae, 2005).In order to characterize multifractal structures within complex nonlinear heterogeneous processes, multifractal detrended fluctuation analysis (MF-DFA; Kantelhardt et al., 2002) was developed based on the framework of DFA, which is mostly used to quantify monofractality.DFA and MF-DFA have been widely applied to evaluate the fractal scaling properties of rainfall and streamflow time series in hydrology (Kantelhardt et al., 2002;Koscielny-Bunde et al., 2006;Labat et al., 2011;Livina et al., 2003;Matsoukas et al., 2000;Zhang et al., 2008).
More specifically, DFA was first adopted in subsurface hydrology by Li and Zhang (2007) to systematically evaluate the fractal dynamics of groundwater systems.They analyzed 4 years of continuous hourly data from seven wells and found that groundwater level fluctuations are likely to follow fractional Brownian motion (fBm) and that temporal scaling crossovers exist in the fluctuations.These findings were later confirmed by Little and Bloomfield (2010), Rakhshandehroo andAmiri (2012), andYu et al. (2016) with the application of DFA to hourly or 15 min interval data for up to 5 years from 7 wells, daily data for 6 years from 2 wells, and daily data from 22 wells that have more than 2500 records, respectively.Rakhshandehroo and Amiri (2012) further utilized MF-DFA to evaluate the multifractality of groundwater level fluctuations and concluded that the extent of multifractality in groundwater level fluctuations is stronger than that in river runoff.
Unlike the general finding of fBm-type behavior in groundwater level fluctuations (Li and Zhang, 2007;Little and Bloomfield, 2010;Rakhshandehroo and Amiri, 2012;Yu et al., 2016), Joelson et al. (2016) found persistent scaling behavior in the analysis of hourly groundwater level fluctuation time series for a 14-month duration and fit the fluctuation data with the Lévy stable distribution to account for the observed non-Gaussian heavy-tailed behavior.
Multiscale multifractal analysis (MMA) was proposed on the basis of MF-DFA, which normally analyzes time series with crossovers only on a predefined large or small scale to obtain the generalized Hurst surface, which simultaneously provides local fractal properties at various scale ranges (Gierałtowski et al., 2012;Wang et al., 2014).To the best of our knowledge, MMA has not yet been applied to analyze time series in hydrology or subsurface hydrology.
In this paper, DFA, MF-DFA, and MMA are applied to systematically evaluate the temporal fractal scaling properties (monofractality and multifractality) of groundwater level fluctuations in two confined aquifer wells with daily data of 70 and 80 years in Texas, USA.Long-term groundwater level data are used, since the Hurst exponent estimated by a larger number of data points tends to be more stable (Weron, 2002).We also check the variation in the estimated Hurst exponent by DFA with different lengths of data and variable time intervals, which is largely unexplored in the aforementioned studies.The possible explanation of the existence of multifractality is studied by MF-DFA and MMA.Furthermore, we investigate the groundwater level fluctuation probability distribution by fitting the data with the α-stable distribution and other distributions, such as Gaussian distribution, gamma distribution, and lognormal distribution, to check if the fBm identified in previous studies is adequate to characterize groundwater level fluctuations.Additionally, we compare the Hurst exponent from fractal analysis with that from the stability index of the fitted α-stable distribution, since the stability index and the Hurst exponent are related under certain conditions (Taqqu et al., 1997).

Methodology
Since the pioneering work of Hurst (1951) on the longmemory behavior (or persistent fractal) of the storage capacity of reservoirs in the Nile River, the Hurst exponent has been regarded as the best-known estimator indicating the magnitude of long-range dependence in time series and has been widely used to study fractal scaling behavior in geophysical sciences, specifically for river flows and turbulence (Nordin et al., 1972;Szolgayova et al., 2014;Vogel et al., 1998), porosity and hydraulic conductivity in subsurface hydrology (Molz and Boman, 1993), climate variability (Bloomfield, 1992;Franzke et al., 2015;Koutsoyiannis, 2003), and sea level fluctuations (Barbosa et al., 2006;Ercan et al., 2013).The Hurst exponent H may be defined as follows: where φ is a given stochastic process, t is time, c is a positive constant, and d is the finite dimension of the time series data; 0 < H < 0.5 demonstrates anti-persistent behavior, H = 0.5 corresponds to uncorrelated noise, 0.5 < H < 1 indicates long-range dependence (i.e., persistent behavior), and H = 1 is for pink noise.
Here, the Hurst exponent was adopted to quantify the scaling properties of groundwater level fluctuation time series.Many methods for the estimation of the Hurst exponent are used in the literature, and different methods may provide significantly different estimates.Detrended fluctuation analysis is chosen here due to its superior performance compared to conventional methods in detecting evolving nonstationarities, which can be very useful to investigate the fractal behavior of time series datasets with different time intervals and in differentiating seasonal trends from the inherent fluctuations of time series (Yu et al., 2016).

Detrended fluctuation analysis
Detrended fluctuation analysis (DFA), also known as variance of the regression residuals, was proposed by Peng et al. (1994).The method is briefly summarized as follows.
First, the original time series {x t } t = 1, 2, . .., n are converted to corresponding sums as (2) Then, {X t } is divided into m (m = n/ l) nonoverlapping blocks Y j of size l, and a least-squares fit (or the local trend) is performed by calculating the variance for each block: where Y l (j ) is the local fitted polynomial trend of first order, second order, or any other higher order.Finally, the root mean-square over all blocks is calculated, yielding the "fluctuation": Fitting a linear line of log(F (l)) against log(l) would indicate the presence of power-law scaling as For fractional Gaussian noise (FGN), α = H , where H is the Hurst exponent.For nonstationary processes (e.g., fractional Brownian motion) α = H + 1 (Heneghan and Mc-Darby, 2000).In this study, the local trend is fitted by a linear line.The DFA method does not assume stationarity in advance.Moreover, it is less sensitive to trends within the data than other approaches, such as the R/S, since a linear regression fit is applied locally in each block.

Multiscale multifractal analysis
Multiscale multifractal analysis (MMA) is a generalization of multifractal detrended fluctuation analysis (MF-DFA), which is developed from DFA (Gierałtowski et al., 2012).In contrast to MF-DFA, which requires the presumption of scaling ranges, MMA is capable of concurrently characterizing different fractal properties (monofractality or multifractality) of time series over a wide range (both small and large) of temporal scales.MMA can be specified as follows.
Based on DFA, the qth-order fluctuation is calculated as (Kantelhardt et al., 2002) If a long-range power-law correlation exists in the time series, then F q (l) for large values of l yields (Kantelhardt et al., 2002) where h(q) is the generalized Hurst exponent and values of h(q) can be interpreted as follows: h ∈ (0, 0.5) indicates antipersistent behavior of the time series, h = 0.5 denotes uncorrelated noise, h ∈ (0.5, 1) indicates persistent behavior of the time series, h = 1.5 corresponds to Brownian motion, and h ≥ 2 indicates black noise; h(q) yields the classical Hurst exponent H when q = 2 for stationary time series and H = h (2) − 1 for nonstationary time series, and h(q) is independent of q for monofractal data and strongly depends on q for time series showing persistent multifractal behavior.The strength of multifractality may be further measured by the Hölder spectrum or singularity spectrum (Feder, 1988).The Hölder exponent α q and the Hölder spectrum (singularity spectrum) f α q can be computed as follows (Kantelhardt et al., 2002): where τ q is the classical multifractal scaling exponent.The strength of multifractality in a time series can be estimated by the width of f α q , which can be illustrated by the range of α q as α q = α max − α min (Koscielny-Bunde et al., 2006).The above estimators show the formulation of MF-DFA.After the calculation of all F q (l) values by MF-DFA, a moving fitting time window, which completely sweeps through the range of scale l along F q (l), is used to study quasicontinuous changes between h(q) dependence and the range of scale l.The fitting procedure is as follows: where f i is a fitting window (i = 1, 2, . .., n) and h f i is the local scaling exponent in f i .For a fixed q, the spectrum of scaling exponents over the whole range of scale l is obtained by h(q, l) = h f 1 , h f 2 , . .., h f n .After plotting the results of h (q, l) for all the q values, the Hurst surface h (q, l), which simultaneously provides the generalized Hurst exponent for multiple scales and q, is obtained (Wang et al., 2014).The capability of MMA, which is inherited from DFA and MF-DFA, is that it can effectively detect observational noise and nonstationarities in time series.Similar to MF-DFA, the results of h (q, l) in MMA characterize large fluctuations in the fragments of data for q > 0, while the results of h (q, l) correspond to small fluctuations for q < 0.

Alpha-stable distributions
The α-stable distributions introduced by Lévy (1925) represent a class of stable laws determined by four parameters: the stability index α, the skewness parameter β, the scale parameter γ , and the location parameter δ.Therefore, the α-stable distribution of a random variable X is usually denoted by X ∼ S α (β, γ , δ).No closed forms exist for the αstable distributions, except for the following three distributions: Gaussian, Cauchy, and Lévy.The α-stable distribution of a random variable, X ∼ S α (β, γ , δ), can be described by the following characteristic function (Samoradnitsky and Taqqu, 1994): where The stability index α is also known as the characteristic exponent and is in the interval of α ∈ (0, 2].The distribution becomes a normal distribution when α = 2.The skewness parameter satisfies −1 ≤ β ≤ 1.The location parameter δ indicates the shift of the peak of the distribution and it is undefined unless α > 1.The distribution is symmetric around δ if β = 0.The scale parameter γ measures the dispersion of the distribution and is always positive (γ > 0).
Stable distributions are heavy tailed, and the tails of these distributions demonstrate asymptotical power law behavior with 0 < α < 2 and −1 < β < 1.One important property of the α-stable distribution is that there is a possible link between the stable distribution and self-affine behavior according to the generalized central limit theorem (Gnedenko and Kolmogorov, 1956).To be more specific, approximation of the tail of the stable distribution X ∼ S α (β, γ , δ) may be shown (Samoradnitsky and Taqqu, 1994): where c α = 1 π sin πα 2 (α).This behavior indicates that αstable distributions can be well accommodated to model selfsimilar processes.The distribution with 1 < α < 2 is of significant interest to researchers as the mean of the distribution can be defined and the variance is infinite.The non-integer α in this range, which is capable of characterizing processes with infinite variance, is related to the Hurst exponent H , presenting long-range dependence and statistical self-similarity properties, as follows (Taqqu et al., 1997): Since the Lévy α-stable distribution is 1/α self-similar, the following equation is also used to describe the relationship between the stability index and the Hurst exponent (Peters, 1994): 3 Data analysis  (more than 1000 days) for Well 1 to become decorrelated, while it takes a couple of years (more than 500 days) for Well 2.Moreover, the ACF plots greatly vary in different 20year intervals of the two datasets (bottom left and right panels in Fig. 2), which may imply that the long-range dependence characteristics of the two wells would vary through time.
The power spectra of Well 1 (1945-2014) and Well 2 (1935-2014) groundwater levels are presented in Fig. 3.The power-law exponents are estimated as 2.44 and 2.08 for Well 1 and Well 2 groundwater levels, respectively, indicating the existence of fractals in both datasets.Hurst exponents can be deduced from the power-law exponents (Heneghan and McDarby, 2000) as 0.72 and 0.54 for Well 1 and Well 2 groundwater levels, respectively.Furthermore, the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test (Kwiatkowski et al., 1992) is conducted to test the stationarity of the data.The null hypothesis for the KPSS test is that a time series is stationary, and the alternative is that data are nonstationary.The estimates of the KPSS statistic are 5.6357 and 1.8012 for Well 1 and Well 2 groundwater levels, respectively, and both reject the null hypothesis at a 1 % significance level, which suggests that the two time series are nonstationary.These results provide reference for the quantification of the Hurst exponent later by DFA.

Monofractal analysis
The Hurst exponents of groundwater level fluctuation data, quantified by the DFA approach over different time intervals, are investigated here.The evolution of the Hurst exponent H through time is shown in Fig. 4, for which the data were chosen in the original order, moving year by year forward from 1945 to 2014 for Well 1 (i.e 1945-1949, 1945-1950, . . . , 1945-2014) and from 1935 to 2014 for Well 2 (i.e 1935-1939, 1935-1940, . . . , 1935-2014). Figure 4a and b    groundwater levels.H is 0.73 and 0.51 for Well 1 and Well 2, respectively, when all the available data are used, which suggests that both datasets indicate long memory.These estimates of Hurst exponents are also consistent with the ones that are deduced from the power-law exponents in Fig. 3.In general, groundwater level fluctuations of Well 1 show persistent fractal behavior (H > 0.5, more specifically H > 0.7) for all investigated time periods, and those of Well 2 vary be-tween persistent and anti-persistent, even showing uncorrelated behavior at certain times.The Hurst exponent H for Well 1 groundwater levels varies between 0.8 and 0.85 for up to 8 years of daily data for end years 1948-1952 (Fig. 4a).It then stabilizes within values of 0.71 and 0.78 for 9 years and longer time durations (for end years after 1952; Fig. 4a).On the other hand, H for Well 2 groundwater levels varies between 0.53 and 0.6 for up to 7 years of daily data for end years 1939-1941 (Fig. 4b).It then stabilizes within 0.46 and 0.53 for time durations longer than 8 years (for end years after 1941; Fig. 4b).As such, the fractal behavior of groundwater levels obtained from shortduration data (in this study, less than 8 years for Well 1 and 7 years for Well 2) may not exhibit stable long-term fractal behavior.These results further imply that the length of time series and the time period covered jointly affect the value of H .The Hurst exponents here demonstrate the ability of DFA in distinguishing the seemingly long-range correlations caused by external effects (such as seasonal trends) from intrinsic fluctuations (Yu et al., 2016), since the ACF plots show very slow decay in both wells (Fig. 2).
Figure 5 presents the Hurst exponents of groundwater level data estimated with different moving time windows (5, 10, and 20 years).Daily data were used in different time windows: 5-year moving window (i.e., 1945-1949, 1946-1950, . . . , 2010-2014), 10-year moving window (i.e., 1945-1954, 1946-1955, . . . , 2005-2014), and 20-year moving window (i.e., 1945-1964, 1946-1965, . . . , 1995-2014). Figure 5a and  b show that the Hurst exponents vary greatly in different time windows (i.e different length of groundwater level fluctuation data) and also do not remain constant even with the same time window when the time window moves in time.Moreover, the results in Fig. 5a and b demonstrate that the Hurst exponent tends to be stable as the time window increases, which is consistent with the results in Fig. 4.
Additionally, the correlation coefficient r is used to investigate the relationship between the Hurst exponent and the variation in groundwater level fluctuations, which is quantified by the coefficient of variation, c v (c v = δ/µ, where δ is the SD of the data and µ is the corresponding average).From Fig. 5c it may be inferred that strong linear correlation exists between c v and H (r > 5), and the correlation becomes stronger as the time window increases from 5 to 20.Meanwhile, for Well 2 groundwater levels the correlation is weaker (r < 0.5), and r increases first and then decreases afterwards following the increase in the time window from 5 to 20 (Fig. 5d).For Well 1 groundwater levels (Fig. 5c), a larger c v normally follows a greater H for 5-, 10-, and 20-year time windows.However, for Well 2 groundwater levels (Fig. 5d), this relationship generally does not hold (especially for the 20-year time window).Figure 5c and d suggest that the variability in groundwater level fluctuations may affect the intrinsic correlation (long memory or short memory) of the data, but it is highly site specific.The different Hurst exponents in different wells may be due to the effect of the heterogeneity of the aquifer materials (Li and Zhang, 2007).
Figure 6 further investigates the variation in the Hurst exponents with box plots for 5-, 10-, 20-, 30-, 40-, and 50-year   .F q as a function of timescale l and the generalized Hurst exponent h as a function of q for the groundwater levels of (a) Well 1 (b) Well 2; (c) the scaling exponent spectrum τ q vs. the moments for the groundwater levels of Well 1 and Well 2; (d) the singularity spectrum for the groundwater levels of Well 1 and Well 2. moving time windows.Unlike the inconsistency of the linear correlation between c v and H , the variation in H in both Well 1 and Well 2 groundwater levels is consistent here.The variation in H for the groundwater levels in both wells generally decreases as the moving time window increases, which confirms the findings in Fig. 5a and b.

Multifractal analysis
The multifractal results obtained by MF-DFA in Fig. 7 include log-log plots of F q (l) against timescale l, the generalized Hurst exponent h(q), the scaling exponent spectrum τ q and the singularity spectrum f α q corresponding to a series of moments q (−5 ≤ q ≤ 5); h(q) is the slope of the linear regression line of the log-log plot for a given q.Clearly, multifractality exists in the groundwater levels of Well 1  and Well 2 (1945-2014), since h(q) greatly varies with q, as demonstrated in Fig. 7a and b, and the relationships between τ q and q in Fig. 6c are not linear for the groundwater levels of Well 1 and Well 2. This also suggests that different exponents should be used to illustrate the fractal scaling behavior (self-affinity) of different time intervals of the data.Moreover, h(q) continuously decreases as q increases in both figures, implying that relatively small fluctuations occur more frequently in the time series than large ones (Grech and Czarnecki, 2009;Rakhshandehroo and Amiri, 2012).
The singularities of the processes in the groundwater levels of Well 1 and Well 2 are revealed in Fig. 7d.The width of the singularity spectrum, α q , is used to measure the level of multifractality.The width of the singularity spectrum α q tends to be zero for monofractal structures and would increase as the level of multifractality of the signal increases; α q was found to be 4.05 for the groundwater levels of Well 1 and 1.07 for the groundwater levels of Well 2. These results indicate a high level of multifractality in both time series, and Well 1 groundwater levels have a stronger multifractality, which further suggests that the multifractal behavior is quite site specific.
Two types of rationale are used to account for multifractality in time series (Kantelhardt et al., 2002).The first type is that a broad probability density function of time series data, which cannot be represented by a regular distribution with finite moments, causes multifractality.The second type is that multifractality is caused by long-range correlations of small and large fluctuations (Kantelhardt et al., 2002;Rakhshandehroo and Amiri, 2012).To distinguish these two types of multifractality, the corresponding randomly shuffled dataset is analyzed.The multifractality will vanish if it is totally due to the second type and will remain otherwise.If the multifractality is due to both types, the shuffled data will present weaker multifractality than the original data (Kantelhardt et al., 2002).
Therefore, a shuffling procedure was conducted to investigate the types of the multifractality for Well 1 and Well 2 groundwater levels.The corresponding multifractality results are shown in Fig. 8.This figure clearly shows that multifractality still exists in the shuffled groundwater level data of Well 1, since dependency between h(q) and q remains (Fig. 8a).The relationship between τ q and q is not linear (Fig. 8c), which further verifies the existence of multifractality in shuffled Well 1 data; α q was 0.18 (Fig. 8d), which indicates a much weaker multifractality compared with α q = 4.05 for the original data.The results for shuffled Well 2 groundwater level data, on the contrary, show that shuffling almost completely destroyed the intrinsic fractal correlations, since h(q) is independent of q (Fig. 8b), τ q is linear with q (Fig. 8c), and the singularity spectrum almost converges to a single point with α q = 0.02 (Fig. 8d), which may indicate an approximate monofractal structure in Well 2 groundwater levels.
The results in Fig. 8 reveal that different types of multifractality exist in Well 1 and Well 2 groundwater level time series.For Well 1, the multifractality is clearly due to the combined effect of a broad probability density function and temporal correlations in diverse magnitudes of fluctuations.Meanwhile, the multifractality is almost purely caused by long-range temporal correlations in small and large fluctuations for Well 2 groundwater levels.
Since the Hurst exponent varies for different time intervals of the groundwater level time series of Well 1 and Well 2 (Figs. 4, 5, and 6) and the small and large fluctuations of temporal correlations contribute to the multifractality of both datasets (Fig. 8), multiscale multifractal analysis (MMA) is adopted to investigate the fractal behavior at different temporal scale ranges in detail, as demonstrated in Fig. 9.It is noted that the generalized Hurst surfaces for the original datasets of both Well 1 and Well 2 groundwater levels (top images in Fig. 9) are far from flat (hill-like shape), which clearly suggests that different fractal scaling exponents are needed to represent fractal behavior at multiple temporal scales for both datasets.In addition, the generalized Hurst exponents at q = 2 are between 1.5 and 2 for Well 1 groundwater levels, indicating persistent behavior, and are mostly within the range between 1 and 1.55 for Well 2 groundwater levels, indicating persistent and anti-persistent fractal behavior (sometimes even uncorrelated).Moreover, the Hurst surfaces for the shuffled time series of Well 1 and Well 2 (bottom images in Fig. 9) show that the surfaces become much flatter than those generated by the original datasets (small fluctuations in the Hurst surfaces still exist after shuffling),which suggests that the shuffling substantially destroys the intrinsic correlations, as consistent with the MF-DFA results.

Relationship between the stability index and the Hurst exponent
Multifractal analysis suggests that multifractality is partially due to a broad probability density distribution that may have infinite moments.However, fBm (fractional Brownian motion), which has been identified as an appropriate model to characterize groundwater level fluctuations (Li and Zhang, 2007;Little and Bloomfield, 2010;Rakhshandehroo and Amiri, 2012;Yu et al., 2016), is Gaussian with finite moments.Therefore, fBm may be inappropriate for the description of physical processes with infinite moments, such as the groundwater level fluctuations in this study.Histograms and normal probability plots for Well 1 and Well 2 groundwater levels in six selected durations of varying length apparently indicate that the Gaussian distribution may not be suitable to represent the groundwater level processes of both wells, especially for Well 1 (Figs. 10 and 11, in which the probability curve would lie on the straight red line if the data are normally distributed).Well 1 groundwater levels clearly show a heavy tail, and Well 2 groundwater levels demonstrate right-skewed behavior.As such, the Lévy alpha-stable distribution, which is non-Gaussian with a heavy tail and has infinite variance, was adopted to fit the groundwater datasets.Moreover, to obtain a relatively comprehensive picture of the underlying probability distribution, the Gaussian distribution, gamma distribution, and lognormal distribution were also used to fit the datasets (the Statistics and Machine Learning Toolbox in Matlab is used for this purpose).The fitting procedure is conducted continuously with the data starting from 1945 for Well 1 and from 1935 for Well 2 and moves forward year by year with the same end year, 2014, for all the fitting durations (at least 15 years of daily data are used to ensure a good characterization of the data distribution).Results for the six selected durations are presented in Figs. 12 and 13 for Well 1 and Well 2 groundwater levels, respectively.
Figure 12 shows that, in general, the Lévy stable distribution fits the groundwater level fluctuation time series of Well 1 over different durations very well.Meanwhile, the www.earth-syst-dynam.net/8/931/2017/Earth Syst.Dynam., 8, 931-949, 2017 other distributions, i.e., normal, gamma, and lognormal, cannot satisfactorily capture the behavior of the groundwater levels of Well 1.This verifies the finding that the irregular distribution of Well 1 groundwater levels contributes to the multifractality.For Well 2 groundwater levels, Gaussian distribution adequately fits the data, except at the peak values (Fig. 13).Furthermore, the fitted stable, gamma, and lognormal distributions converge to the Gaussian distribution.This may imply that fBm may partially represent the behavior of Well 2 groundwater levels, which has the Hurst exponent fluctuating between 0.48 and 0.52 (Fig. 14b).Furthermore, the stability index α of the stable distribution is related to the Hurst exponent H given by a relationship between α and H . Two formulae (Eqs.12 and 13) are used to estimate H .The estimated H is then compared with that estimated by DFA (Fig. 14).With respect to the difference between the Hurst exponent estimated by DFA (for both Well 1 and Well 2 groundwater levels) and that deducted from either H = 1 α or H = 3−α 2 , the relative difference is less than 10 % (even less than 1 % in some time intervals) for most of the comparisons.For Well 1 groundwater levels, the Hurst exponent by H = 1 α generally matches better with H estimated by DFA than that by H = 3−α 2 (Fig. 14a), although for some durations, such as 1950-2014 or 1955-2014, the latter performs better than the first.Figure 14c further shows that the stability index α is strongly correlated with the coefficient of variation c v of the groundwater level fluctuation data from Well 1, since the correlation coefficient is as high as −0.84.It suggests that a larger c v of Well 1 groundwater levels would probably imply a smaller α.A stability index α = 2 for all the stable distributions for the groundwater levels data from Well 2 (α = 2 corresponds to Gaussian distribution) is found.This may be due to the fact that the Lévy stable distribution here is restricted in the range 1 < α ≤ 2, which corresponds to 0.5 ≤ H < 1.However, Well 2 groundwater levels do not have long memory in some time intervals.The relationship between H and α would completely fail when H < 0.5.However, the resulting stability index α = 2 for Well 2 groundwater levels is acceptable considering that the difference between H estimated by DFA and H estimated by the stability index is less than 5 % for all the time periods.This result is also consistent with Fig. 11 in which Gaussian distribution is capable of capturing the main groundwater level fluctuation patterns of Well 2.
The results indicate that fBm, which has Gaussian characteristics, may be a reasonable model for representing groundwater level fluctuations under certain conditions, such as in the case of the dataset of Well 2, which has the Hurst exponent fluctuating close to around 0.5.However, fBm may be an insufficient model for capturing the behavior of groundwater fluctuations in other cases, for example in the case of the groundwater levels of Well 1, in which a non-Gaussian distribution, such as a heavy-tailed stable distribution (Lévy motion), is needed instead.In the presence of long memory, fractional Lévy motion may be more appropriate to model and forecast the groundwater dynamics.
It is important to note that the results obtained so far are limited to the analysis of temporal correlations of the groundwater level fluctuations at certain locations.The properties of the groundwater levels at the two wells, such as their fractal behavior and underlying distributions, are highly different from each other, which confirms that the results are site specific.Well 1 and Well 2 are chosen because the groundwater level fluctuation records of these two wells are long and complete.In addition, these two wells are very representative in terms of fractal scaling behavior and underlying probability density distribution.
Groundwater dynamics in aquifers result from multiple complex dynamic processes, such as hydrologic processes (precipitation, river runoff, etc.), the hydraulic properties of soil and aquifers, and anthropogenic perturbations (such as the construction of reservoirs and pumping of water).These processes and properties vary at different spatiotemporal scales, which directly or indirectly affect groundwater systems.The two confined aquifer wells analyzed in this study are located at the same type of aquifer but present drastically different dynamics of groundwater level fluctuations.The results obtained herein may be attributed to the time-space heterogeneity of aquifer characteristics, hydrometeorological conditions, and even anthropogenic forcing, but detailed research, such as the employment of time-space analysis, needs to be conducted to justify this and to account for the effect of heterogeneity on fractal behavior at different temporal scales.The non-Gaussian fractal property of the groundwater system in Well 1 that demonstrates long memory provides further insight for the resulting transport processes in the porous medium, which may also present non-Gaussian features with memory similar to the non-Gaussian behavior that is found in the precipitation time series in other studies (Joelson et al., 2016;Lovejoy and Mandelbrot, 1985).Unlike Well 1 groundwater levels, the origin of multifractality for the Well 2 groundwater levels is difficult to explain due to the very weak multifractality after the shuffling.An intuitive explanation may be that it is due to noise.However, the fractal structure is not affected by dynamical noise (Serletis, 2008).Additionally, Gaussian distribution may partially represent the dataset of Well 2 groundwater levels, but it fails to capture the peak of the skewed distribution of Well 2 groundwater levels, which may imply that an irregular distribution that also holds certain Gaussian characteristics may be needed to fully characterize the groundwater dynamics of Well 2.

Conclusions
In this study, the fractal scaling properties of groundwater level fluctuations of two confined aquifer wells with 70 and 80 years of daily data were analyzed.Detrended fluctuation analysis (DFA) was utilized to quantify the Hurst exponent and monofractality.The DFA results indicated that fractals exist in the groundwater level time series of both wells, and it was shown that the Hurst exponent is highly dependent on the length and specific time period of the time series.A persistent scaling pattern was found for all investigated time periods of Well 1 groundwater levels (Hurst exponent, H > 0.5), and the scaling pattern varied between anti-persistent and persistent regimes for Well 2 groundwater levels.The Hurst exponent H for Well 1 groundwater levels fluctuated between 0.8 and 0.85 for up to 8 years of daily data for end years 1948-1952 (Fig. 4a) and then stabilized within the range of 0.71-0.78for 9 years and longer time durations (for end years after 1952; Fig. 4a).On the other hand, H for Well 2 groundwater levels fluctuated between 0.53 and 0.6 for up to 7 years of daily data for end years 1939-1941 (Fig. 4b) and then stabilized within the range 0.46-0.53for durations longer than 8 years (for end years after 1941; Fig. 4b).
Multifractal detrended fluctuation analysis (MF-DFA) and multiscale multifractal analysis (MMA) were adopted to examine the multifractality and multifractal behavior at different temporal scales for confined groundwater levels.Although the MF-DFA results showed that a relatively high level of multifractality exists for the groundwater levels of both wells, a stronger multifractality was observed for the dataset of Well 1 compared to that of Well 2. The observed multifractality was postulated to originate from the combined effect of the underlying irregular probability distributions and different magnitudes of fluctuations in multiple longrange temporal correlations for Well 1 groundwater levels and mostly long-range temporal correlations in small and large fluctuations for Well 2 groundwater levels.Moreover, the MMA results confirmed the existence of multifractality and diverse correlations of groundwater levels over different timescales.For Well 1, the  generally matches H estimated by DFA than better that by H = 3−α 2 .The stability index α is strongly correlated with the coefficient of variation c v of the groundwater level fluctuation data from Well 1, since the correlation coefficient is as high as −0.84.For Well 2, a stability index α = 2 for all the fitted stable distributions for the groundwater level data from Well 2 (α = 2 corresponds to Gaussian distribution) is found.
Furthermore, the underlying probability distribution of groundwater level fluctuations for Well 1 represented mainly long-memory characteristics, which were fitted reasonably well by the Lévy stable distributions for various time periods.On the other hand, those of Well 2 represented mainly Gaussian characteristics, which sometimes failed to capture the peaks of the skewedprobability distributions of Well 2 groundwater levels.Time series analysis of groundwater level fluctuations of the two wells demonstrated that the observed fractal behavior is site specific, and there is a need for generalized governing equations of groundwater flow processes that can model both the long-memory behavior and the Brownian finite-memory behavior (Kavvas et al., 2017;Tu et al., 2017).

Appendix A: Selection criteria for the two wells
As of May 2016, there are 257 monitoring wells, both active and inactive, reported on the web page of Water Data for Texas (http://waterdatafortexas.org/groundwater/).The groundwater level datasets consist of data taken from 77 confined aquifer wells, 179 unconfined aquifer wells, and 1 unknown aquifer well.Since the focus of this study is on the confined aquifer well, the spatial distribution of the 77 datasets related to the confined aquifer is provided in Fig. A1.
The longest dataset has more than 81 years of record with an approximate 2.6 % missing rate, and the second-longest dataset includes more than 72 years of record with an approximate 3.6 % missing rate.The third-longest dataset has more than 10 % missing measurements and has about one-third of the length of the second-longest dataset.Therefore, the first-and second-longest groundwater level records were analyzed in this study.Record lengths (in days) and percentage of missing data for the 20 longest groundwater level records reported by Water Data for Texas are presented in Fig. A2.
It can be seen from Fig. A2 that only the first two records are of excellent data quality with respect to length and completeness.Therefore, the groundwater level fluctuations of these two confined aquifer wells are analyzed in this study.The results indicate two different behaviors of the groundwater level fluctuations, i.e., Gaussian and non-Gaussian, which are not reported or compared in previous studies on the fractal scaling analysis of groundwater level fluctuations.Therefore, the results of this behavior with respect to these two confined aquifer wells are reported in the paper.These two wells are both located at the Edwards (Balcones Fault Zone) aquifer, which primarily consists of partially dissolved limestone.However, the dynamics of the groundwater level fluctuations in these two wells behave drastically differently, which may imply rather different climate-related and anthropogenic perturbations in these two wells.Unfortunately, due to the lack of high-quality datasets and detailed information about the aquifers in this area, further discussion on the regionalization of the fractal properties is difficult.

Figure 1 .
Figure 1.Groundwater level time series data of (a) Well 1 and (b) Well 2.
clearly show that the Hurst exponent varies through time for both Well 1 and Well 2 groundwater levels.Box plots in Fig. 4c demonstrate that the mean and variance of the Hurst exponent through time differ noticeably for both Well 1 and Well 2 www.earth-syst-dynam.net/8/931/2017/Earth Syst.Dynam., 8, 931-949, 2017

Figure 4 .
Figure 4. (a) The evolution of the Hurst exponent through time for Well 1 groundwater levels, which starts from 1945; (b) the evolution of the Hurst exponent through time for Well 2 groundwater levels, which starts from 1935; (c) box plots of the Hurst exponents in panels (a) and (b).

Figure 5 .
Figure 5. (a) The Hurst exponent of Well 1 groundwater levels estimated by DFA within different time windows.(b) The Hurst exponent of Well 2 groundwater levels estimated by DFA within different time windows.(c) The relationship between the coefficient of variation and the Hurst exponent obtained in panel (a).(d) The relationship between the coefficient of variation and the Hurst exponent obtained in panel (b).

Figure 6 .
Figure 6.Box plots of the Hurst exponents under different moving time windows for (a) Well 1 groundwater levels and (b) Well 2 groundwater levels.

Figure 7
Figure7.F q as a function of timescale l and the generalized Hurst exponent h as a function of q for the groundwater levels of (a) Well 1 (b) Well 2; (c) the scaling exponent spectrum τ q vs. the moments for the groundwater levels of Well 1 and Well 2; (d) the singularity spectrum for the groundwater levels of Well 1 and Well 2.

Figure 8 .
Figure8.F q as a function of timescale l and the generalized Hurst exponent h as a function of q after shuffling for (a) Well 1 groundwater levels and (b) Well 2 groundwater levels.(c) The scaling exponent spectrum τ q vs. the moments for Well 1 and Well 2 groundwater levels after shuffling and (d) the singularity spectrum for Well 1 and Well 2 groundwater levels after shuffling.

Figure 9 .
Figure 9. Multiscale multifractal analysis for (a) Well 1 groundwater levels and (b) Well 2 groundwater levels.The top figure (in red, yellow, and green) is for the original data and the bottom (in blue) is for the shuffled data.The thick black line indicates MF-DFA results for a given temporal scale, and the solid dots show the generalized Hurst exponents at q = 2 by DFA over different scales.

Figure 10 .
Figure 10.Histograms and normal probability plots for various time intervals of groundwater levels of Well 1.

Figure 11 .
Figure 11.Histograms and normal probability plots for various time intervals of groundwater levels of Well 2.

Figure 12 .
Figure 12.Probability density function and cumulative distribution function (first two rows and last two rows, respectively) of groundwater level fluctuation time series data of Well 1.

Figure 13 .
Figure 13.Probability density function and cumulative distribution function (first two rows and last two rows, respectively) of groundwater level fluctuation time series data of Well 2.

Figure 14 .
Figure 14.(a) Well 1 groundwater levels: comparison between the Hurst exponent estimated by DFA and by stability index.(b) Well 2 groundwater levels: comparison between the Hurst exponent estimated by DFA and by stability index.(c) Well 1 groundwater levels: the coefficient of variation c v vs. the stability index obtained in panel (a).(d) Well 2 groundwater levels: the coefficient of variation c v vs. the stability index obtained in panel (b).

Figure A1 .
Figure A1.Spatial distribution of the confined aquifer wells in Texas, USA reported by Water Data for Texas.The two red stars denote the locations of the wells that have the first-and secondlongest records; yellow solid circles denote the locations of the other 75 wells.

Figure A2 .
Figure A2.Record lengths (in days) and percentage of missing data for the 20 longest groundwater level records reported in Texas by Water Data for Texas.

Table 1 .
Statistics and geophysical properties of studied wells in Texas, USA.
The autocorrelation function (ACF) in Fig.2shows very slow decay in both datasets, and the dataset of Well 1 decays more slowly than that of Well 2. In fact, it takes several years Earth Syst.Dynam., 8, 931-949, 2017www.earth-syst-dynam.net/8/931/2017/Note: mean represents the mean groundwater level (hydraulic head) depth below land surface.