Evaluation of statistical methods for quantifying fractal scaling in water-quality time series with irregular sampling

River water-quality time series often exhibit fractal scaling, which here refers to autocorrelation that decays as a power law over some range of scales. Fractal scaling presents challenges to the identification of deterministic trends because (1) fractal scaling has the potential to lead to false inference about the statistical significance of trends and (2) the abundance of irregularly spaced data in water-quality monitoring networks complicates efforts to quantify fractal scaling. Traditional methods for estimating fractal scaling – in the form of spectral slope (β) or other equivalent scaling parameters (e.g., Hurst exponent) – are generally inapplicable to irregularly sampled data. Here we consider two types of estimation approaches for irregularly sampled data and evaluate their performance using synthetic time series. These time series were generated such that (1) they exhibit a wide range of prescribed fractal scaling behaviors, ranging from white noise (β = 0) to Brown noise (β = 2) and (2) their sampling gap intervals mimic the sampling irregularity (as quantified by both the skewness and mean of gap-interval lengths) in real water-quality data. The results suggest that none of the existing methods fully account for the effects of sampling irregularity on β estimation. First, the results illustrate the danger of using interpolation for gap filling when examining autocorrelation, as the interpolation methods consistently underestimate or overestimate β under a wide range of prescribed β values and gap distributions. Second, the widely used Lomb–Scargle spectral method also consistently underestimates β. A previously published modified form, using only the lowest 5 % of the frequencies for spectral slope estimation, has very poor precision, although the overall bias is small. Third, a recent wavelet-based method, coupled with an aliasing filter, generally has the smallest bias and root-meansquared error among all methods for a wide range of prescribed β values and gap distributions. The aliasing method, however, does not itself account for sampling irregularity, and this introduces some bias in the result. Nonetheless, the wavelet method is recommended for estimating β in irregular time series until improved methods are developed. Finally, all methods’ performances depend strongly on the sampling irregularity, highlighting that the accuracy and precision of each method are data specific. Accurately quantifying the strength of fractal scaling in irregular water-quality time series remains an unresolved challenge for the hydrologic community and for other disciplines that must grapple with irregular sampling.

where C1 is a constant, is satisfied by some ∈ (0,1) (Boutahar et al., 2007;Beran et al., 2013). 70 Equivalently, Xt is a long-memory process if, in the spectral domain, the condition 71 →0 | | ( ) = 2 > 0 (3) is satisfied by some ∈ (0,1), where C2 is a constant and ( ) is the spectral density function 72 of Xt, which is related to ACF as follows (which is also known as the Wiener-Khinchin theorem): where is angular frequency (Boutahar et al., 2007). 74 One popular model for describing long-memory processes is the so-called fractional auto-75 regressive integrated moving-average model, or ARFIMA (p, q, d), which is an extension of 76 ARMA models and is defined as follows:

82
In addition to a slowly decaying ACF, a long-memory process manifests itself in two other 83 equivalent fashions. One is the so-called "Hurst effect", which states that, on a log-log scale, the 84 range of variability of a process changes linearly with the length of the time period under 85 consideration. This power-law slope is often referred to as the "Hurst exponent" or "Hurst 86 coefficient" H (Hurst, 1951), which is related to d as follows: spectra that are linear on log-log axes (Lomb, 1976;Scargle, 1982;Kirchner, 2005). 91 Mathematically, this inverse proportionality can be expressed as: where 3 is a constant and the scaling exponent β is termed the "spectral slope." In particular, for 6 methods have been extensively adopted, they are unfortunately only applicable to regular (i.e., 116 evenly spaced) data, e.g., daily streamflow discharge, monthly temperature, etc. In practice, 117 many types of hydrological data, including river water-quality data, are often sampled irregularly 118 or have missing values, and hence their strengths of fractal scaling cannot be readily estimated 119 with the above traditional estimation methods. 120 Thus, estimation of fractal scaling in irregularly sampled data is an important challenge for 121 hydrologists and practitioners. Many data analysts may be tempted to interpolate the time series 122 to make it regular and hence analyzable (Graham, 2009). Although technically convenient, 123 interpolation can be problematic if it distorts the series' autocorrelation structure (Kirchner and 124 Weil, 1998). In this regard, it is important to evaluate various types of interpolation methods 125 using carefully designed benchmark tests and to identify the scenarios under which the 126 interpolated data can yield reliable (or, alternatively, biased) estimates of spectral slope. 127 Moreover, quantification of fractal scaling in real-world water-quality data is subject to 128 several common complexities. First, water-quality data are rarely normally distributed; instead, 129 they are typically characterized by log-normal or other skewed distributions (Hirsch et al., 1991;130 Helsel and Hirsch, 2002), with potential consequences for β estimation. Moreover, water-quality 131 data also tend to exhibit long-term trends, seasonality, and flow-dependence (Hirsch et al., 1991;132 Helsel and Hirsch, 2002), which can also affect the accuracy of β estimates. Thus, it may be 133 more plausible to quantify β in transformed time series after accounting for the seasonal patterns 134 and discharge-driven variations in the original time series, which is also the approach taken in 135 this workpaper. For the trend aspect, however, it remains a puzzle whether the data set should be 136 de-trended before conducting β estimation. Such de-trending treatment can certainly affect the 137 estimated value of β and hence the validity of (or confidence in) any inference made regarding 138 the statistical significance of temporal trends in the time series. This somewhat circular issue is 139 beyond the scope of our current work --it has been previously discussed in the context of short-140 term memory (Zetterqvist, 1991 may be broadly applicable to irregularly sampled data in many other scientific disciplines.

168
The rest of the paper is organized as follows. We propose a general approach for modeling 169 sampling irregularity in typical river water-quality data and discuss our approach for simulating 170 irregularly sampled data (Section 2). We then introduce the various methods for estimating 171 fractal scaling in irregular time series and compare their estimation performance (Section 3). We 172 close with a discussion of the results and implications (Section 4).

175
River water-quality data are often sampled irregularly. In some cases, samples are taken 176 more frequently during particular periods of interest, such as high flows or drought periods; here 177 we will address the implications of the irregularity, but not the (intentional) bias, inherent in such 178 a sampling strategy. In other cases, the sampling is planned with a fixed sampling interval (e.g.,

219
To visually illustrate these gap distributions, representative samples of irregular time series 220 are presented in Figure 1 for the three special processes described above (Section 1.2), i.e.,

230
The above modeling approach was applied to real water-quality data from two large river variance of Δt * as follows: Alternatively, a maximum likelihood approach can be used, which employs the "fitdist" function 242 in the "fitdistrplus" R package (Delignette-Muller and Dutang, 2015). In general, the two 243 approaches have produced similar results, which are summarized in Table 1

255
For the Lake Erie and Ohio tributary monitoring program (6 sites), the records of nitrate-256 plus-nitrite (NOx) and TP were examined. According to the maximum likelihood approach, the 257 shape parameter λ is approximately 0.01 for both constituents (

264
To evaluate the various β estimation methods, our first step was to use Monte Carlo 265 simulation to produce time series that mimic the sampling irregularity observed in real water-266 quality monitoring data. We began by simulating regular (gap free) time series using the

322
B4 and B5 were implemented using the "na.locf" function in the "zoo" R package (Zeileis and package (R Development Core Team, 2014). An illustration of these interpolation methods is 325 provided in Figure 4. The interpolated data, along with the original regular data (designated as 326 A1) were analyzed using the Whittle's maximum likelihood method for β estimation, which was 327 implemented using the "FDWhittle" function in the "fractal" R package (Constantine and The second type of approaches estimates β directly from in the irregularly sampled data estimated and the spectral slope is fit using all frequencies (Lomb, 1976;Scargle, 1982).

335
This is a classic method for examining periodicity in irregularly sampled data, which is weighted wavelet spectrum (Foster, 1996) to suppress spectral leakage from low 348 frequencies and applies an aliasing filter (Kirchner, 2005) to remove spectral aliasing 349 artifacts at high frequencies.

350
C1a was implemented using the "spec.ls" function in the "cts" R package (Wang, 2013). C2 was 351 run in C, using codes modified from those in Kirchner and Neal (2013). should lead to spectral slopes that are flatter than expected (Kirchner, 2005), and thus to 367 underestimates of LRD.

368
For the simulated irregular data, the estimation methods differ widely in their performance.  and precision when λ is 1 or 10, a slightly better performance when λ is 0.1, but a worse 393 performance when λ is 0.01.

394
The shape parameter λ greatly affects the performance of the estimation methods. All the  However, C1b overall shows less severe bias than C1a. Finally, the wavelet method (C2) shows 441 generally the smallest bias among all methods. However, its performance advantage is not as great when the time series has small λ values (i.e., very skewed gap intervals), as noted above, 443 which may be due to the fact that the aliasing filter was designed for regular time series. In terms 444 of SD (Figure S7), method C1b performs the worst among all methods (as noted above), method 445 B6 and B8 perform poorly for large βtrue values, and method C2 performs poorly for βtrue = 0. In 446 terms of RMSE (Figure 7), methods B1, B8, C1a, and C1b perform well for small βtrue values

461
In this section, the proposed estimation approaches were applied to quantify β in real water-462 quality data from the two monitoring programs presented in Section 2.2 ( Table 1)  The estimated β values for the concentration residuals are summarized in Figure 10. Clearly, 489 the estimated β varies considerably with the estimation method. In addition, the estimated β 490 varies with site and constituent (i.e., TP, TN, or NOx.) Our discussion below focuses on the 491 wavelet method (C2), because it is established above that this method performs better than the 492 other estimation methods under a wide range of gap conditions. We emphasize that it is beyond 493 our current scope to precisely quantify β in these water-quality data sets, but our simulation 494 results presented above (Section 3.2) can be used as references to qualitatively evaluate the 495 reliability of C2 and/or other methods for these data sets.

496
For TN and TP concentration data at the Chesapeake River Input Monitoring sites (Table 1), 497 μ varies between 9.5 and 24.4, whereas λ is ~1.0. Thus, the simulated gap scenario of NB(μ = 14, 498 λ = 1) can be used as a reasonable reference to assess methods' reliability (Figure 8). Based on 499 method C2, the estimated β ranges between β = 0.36 and β = 0.61 for TN and between β = 0.30 500 and β = 0.58 for TP at these sites (Figure 10). For such ranges, the simulation results indicate that method C2 tends to moderately under-estimate β under this gap scenario (Figure 8), and 502 hence spectral slopes for TN and TP at these Chesapeake sites are likely probably slightly higher 503 than those presented above.

504
For NOx and TP concentration data at the Lake Erie and Ohio sites (Table 1)

531
River water-quality time series often exhibit fractal scaling behavior, which presents highlighting that the accuracy and precision of each method are data-specific.

552
Overall, these results provide new contributions in terms of better understanding and 553 quantification of the proposed methods' performances for estimating the strength of fractal 554 scaling in irregularly sampled water-quality data. In addition, the work has provided an 555 innovative and general approach for modeling sampling irregularity in water-quality records.

556
Moreover, this work has proposed and demonstrated a generalizable framework for data 557 simulation (with gaps) and β estimation, which can be readily applied toward the evaluation of to 558 evaluate other methods that are not covered in this work. More generally, the findings and 559 approaches may also be broadly applicable to irregularly sampled data in other scientific 560 disciplines. Last but not least, we note that accurate quantification of fractal scaling in irregular water-quality time series remains an unresolved challenge for the hydrologic community and for 562 many other disciplines that must grapple with irregular sampling.  processes that correspond to white noise (β = 0), pink noise (β = 1), and Brown noise (β = 2).

707
The 1 st row shows the simulated time series without any gap. The three rows below show the 708 same time series as in the 1 st row but with data gaps that were simulated using three different 709 negative binomial (NB) distributions, that is, 2 nd row: NB(λ = 1, μ = 1); 3 rd row: NB(λ = 1, μ =  Table 1 720 for estimated mean and shape parameters.    Table 1 for site and data details.