Ubiquitous Fractal Scaling and Filtering Behavior of Hydrologic Fluxes and Storages from A Mountain Headwater Catchment

Abstract

We used the weighted wavelet method to perform spectral analysis of observed long-term precipitation, streamflow, actual evapotranspiration, and soil water storage at a sub-humid mountain catchment near Tucson, Arizona, USA. Fractal scaling in precipitation and the daily change in soil water storage occurred up to a period of 14 days and corresponded to the typical duration of relatively wet and dry intervals. In contrast, fractal scaling could be observed up to a period of 0.5 years in streamflow and actual evapotranspiration. By considering long-term observations of hydrologic fluxes and storages, we show that, in contrast to previous findings, the phase relationships between water balance components changed with component period and were not perfectly in or out of phase at all periods. Self-averaging behavior was apparent, but the temporal scales over which this behavior was applicable differed among the various water balance components. Conservative tracer analysis showed that this catchment acted as a fractal filter by transforming white noise in the precipitation input signal to a 1/f flicker in the streamflow output signal by means of both spatial and temporal subsurface advection and dispersion processes and soil wetting properties. This study provides an improved understanding of hydrological filtering behavior in mountain critical zones that are critical sources of water and ecosystem services throughout the world.

1. Introduction

Fractal scaling behavior has been broadly observed in environmental research [1,2,3] and across scientific disciplines (e.g., [4,5,6]). In hydrology, catchment-scale water balance studies have found widespread evidence of fractal scaling in the relationship between precipitation and streamflow as input and output signals, respectively [7,8,9,10]. However, less is known about whether such scaling behavior is generally applicable to all components of the catchment water balance, and if so, at what temporal scale(s) fractal scaling is likely to apply.

Fractal scaling behavior provides information about relationships among water balance components. Previous studies have investigated this type of scaling using global circulation or land-atmosphere-ocean system models to simulate river basin discharge [8,9,10,11]. In contrast, Kirchner, et al. [7] characterized fractal scaling behavior of streamflow and precipitation fluxes using high-density long-term observations for humid catchments at Plynlimon, Wales, U.K., where fractal behavior was shown to occur for periods less than ~1 year, but not at longer periods. Notwithstanding, it remains unclear whether the same fractal scaling behavior can be expected for sites differing in climate and geology from the Plynlimon catchments, and if there is evidence for similar fractal scaling behavior in other components of the water balance.

Time-series data form the basis of research into the dynamic relationships between various water balance components. For example, Bales, et al. [12] suggested that evapotranspiration losses are local with respect to precipitation, i.e., that an in-phase relationship exists between the two components, and McDonnell [13] showed that evapotranspiration water losses are out of phase with precipitation (see also [14]), both in California, USA. However, hydrological fluxes from a catchment (evapotranspiration and streamflow) are expected to sample subsurface flow paths of different lengths and transit times [15,16], leading to dynamic, rather than constant, phase relationships. As a result, there is a need for research addressing the dynamic phase relationship hypothesis between various water balance components.

Hydrologists traditionally assume that the long-term average behavior of hydrologic fluxes and/or storages are stable estimates of those processes. In this way, they implicitly assume self-averaging type behavior, and therefore use the long-term averages to characterize catchment behavior (e.g., Sanford, et al. [17]). However, Kirchner, et al. [7,18], and Kirchner and Neal [19] have shown that chemical signatures in stream waters can be fractal in nature, i.e., that there is no characteristic time-scale that explains the implied variability in tracer concentration data. Accordingly, a growing body of literature suggests the potential for non-self-averaging type behavior of chemical constituents in stream water. It is thus logical to investigate whether the non-self-averaging behavior of a conservative tracer is caused by non-self-averaging behavior of the subsurface storages and hydrologic fluxes that are responsible for storage and discharge of the tracer. This question may be addressed by examining the occurrence of temporal self-averaging behavior across water balance components.

Godsey, et al. [20] reported a one-over-frequency (1/f) signal in concentration time series data of conservative solutes in stream water from humid and mostly coastal catchments with maximum elevations of ~580 m above sea level (asl). Kirchner and Neal [19] also reported fractal scaling behavior with a 1/f signal for several reactive and non-reactive solutes in the humid, low elevation Plynlimon, U.K., catchments; however, it has not been established whether fractal scaling behavior is also applicable to higher elevation mountain catchments. Kirchner and Neal [19] suggested that the 1/f signal at Plynlimon could be explained by catchment-scale advection and dispersion processes. Alternatively, Harman [21] indicated that catchment-scale transient storage effects may also explain 1/f signals in the same time series results. In either case, more research is needed to evaluate the potential causes of temporal transformations in non-reactive solutes from precipitation to streamflow.

In order to improve understanding of fractal scaling behavior at the catchment scale, this study used long-term daily observations of hydrologic fluxes and storages from a headwater mountain catchment in Arizona, USA, to address the following questions:

(1) Is fractal behavior pervasive across all components of the catchment-scale water balance and across all time periods?

(2) Do phase relationships between various water balance components suggest period-invariant behavior?

(3) Do long-term observations of various water balance components support temporal self-averaging behavior?

(4) What processes control the input versus output concentrations of a conservative tracer?

Accordingly, this research applied power spectral and root mean square error analyses to long-term daily observations of the water balance including cross-spectral analysis between pairs of components. This analysis demonstrates ubiquitous fractal scaling behavior up to a certain period after which point it breaks down. It also suggests dynamic phase relationships between the various water balance components that can evolve with component period. This work contributes to an improved understanding of hydrological fractal scaling behavior in mountainous headwater catchments that are disproportionally important to the sustainability of global water sources [22,23,24,25].

2. Study Site and Data

2.1. Study Site

The study site, Marshall Gulch catchment (MGC), is located within the Santa Catalina Mountains Critical Zone Observatory near Tucson, Arizona, USA (32°25′38.46″ N, 110°45′28.28″ W; Figure 1A). The surface elevation of MGC varies from ~2285 m above sea level (asl) to ~2632 m asl, with a mean elevation of 2428 m asl and a mean topographic gradient of ~22 degrees. MGC is mostly composed of granite and micaceous schist at upper and lower elevations, respectively [26]. The soil is mostly sandy loam [27] with soil depth varying from 0 to 1.5 m [28]. Soils overlying schist are thicker with greater clay content than soils overlying granite [29]. Vegetation at MGC is predominantly composed of Rocky Mountain aspen forest and woodland and Madrean upper montane conifer-oak forest and woodland at upper elevations, and Madrean pine-oak forest and woodland at lower elevations (GIS data source: [30]). Note that ~64% of the site was affected by fire (Aspen fire) in 2003 [31]. Between 1981 and 2010, the mean annual precipitation and air temperature at MGC were 920 mm and 10.8 °C [32]; mean annual precipitation during the study period (2008 to 2017) was 654 mm [15]. Shallow groundwater in MGC follows multiple flow paths through soil and fractured bedrock, as shown schematically in Figure 1B. Flow path activity depends on seasonal conditions and the state of shallow storage [16].

2.2. Long-Term Observations of Catchment-Scale Hydrologic Fluxes and Soil Water Storage

We used ten water years (WY; water year n is defined as 1 July of calendar year n − 1 to 30 June of year n) of catchment-scale daily hydrologic fluxes and soil water storage data (Figure 2A) to examine potential fractal scaling behavior of the water balance between 2008 and 2017. For this analysis, we specifically considered precipitation, streamflow, actual evapotranspiration (AET), and soil water storage. Unless otherwise stated, information presented in this section is based on Dwivedi, et al. [15].

2.2.1. Precipitation (P)

Seven rain gages- Double Rock (2338 m asl), Fern Valley (2381 m asl), Schist (2338 m asl), Outlet (2325 m asl), Burnt Meadow (2396 m asl), Granite (2378 m asl) and MG-Weir (2283 m asl)—measure spatial variability in rainfall across MGC at 15 minute temporal resolution (Figure 1A). Each gage has a cluster of three tipping bucket range gages (Rainew 111) except for the MG-Weir gage, which has a single tipping bucket. Because these seven gages do not cover the highest elevations at MGC, an additional precipitation time series from the Mt. Lemmon rain gage was used to represent high-elevation precipitation (Figure 1A). The Mt. Lemmon rain gage records precipitation at 15 minute intervals and is a heated gage (Campbell Scientific 385). The precipitation time series at each station was transformed into daily precipitation by summing any consecutive ninety-six 15 minute time steps. The daily precipitation records at all eight stations were combined into a catchment-scale daily precipitation time series using Thiessen polygon-derived weighting (Figure 2A) [15].

2.2.2. Streamflow (Q)

Streamflow at MGC was measured at 30 minute intervals at the MG-Weir site (Figure 1A) using a pressure transducer (U20-001-01; Onset) with a maximum error of 0.62 kPa and an accuracy of 0.02 kPa, and a previously derived stagedischarge relationship [34]. Measured streamflow (L/s) was converted to mm/30 minute and averaged to a daily time step (Figure 2B).

2.2.3. Actual Evapotranspiration (AET)

Dwivedi, et al. [15] estimated daily AET values using tower-based eddy covariance data from the Mt. Bigelow (US-MtB) AmeriFlux site, located ~1.5 km east of MGC at 2573 m asl (Figure 1) [35]. Because the potential exists for a spatial scale mismatch between the MGC and the flux tower footprint, and the vegetation on Mt. Bigelow was not affected by the Aspen fire, the flux tower evapotranspiration time series data were corrected using a water balance method (similar to [36]) by applying Equation (1) below at the catchment-scale.

$$\left(\mathrm{P}\left(\mathrm{t}\right)-\mathrm{Q}\left(\mathrm{t}\right)-\mathrm{AET}\left(\mathrm{t}\right)\right)\mathsf{\Delta }\mathrm{t}-{\mathrm{S}}_{\mathrm{soil}}\left(\mathrm{t}\right)-{\mathrm{S}}_{\mathrm{FBR}}\left(\mathrm{t}\right)=0,$$

(1)

where Δt is set to 1 day, P(t), Q(t), AET (t), $\mathsf{\Delta }{\mathrm{S}}_{\mathrm{soil}}\left(\mathrm{t}\right)$ and $\mathsf{\Delta }{\mathrm{S}}_{\mathrm{FBR}}\left(\mathrm{t}\right)$ are catchment-scale precipitation, streamflow, actual ET, and changes in soil and fractured bedrock water storages, respectively, and it is assumed that the total change in fractured bedrock water storage is zero for any water year. Due to the availability of the observed precipitation, streamflow and soil water storage data, the daily catchment-scale AET data were also estimated for water year 2008 at MGC (Figure 2C) [15].

2.2.4. Soil Water Storage (S_soil)

There are eight soil moisture sensor locations within the Schist sub-catchment and three soil moisture sensor locations within the Granite sub-catchment (Figure 1A). At each sensor location, soil moisture was measured using five vertically-distributed sensors (EC-5; Decagon Devices Inc.). The sensor resolution is 0.001 m³/m³ in mineral soils and the sensor accuracy is ±0.03 m³/m³ using the factory calibration for mineral soils [37]. No sensors were replaced during the study period in order to avoid disturbance to nearby sensors. The volumetric soil moisture data were recorded at 30 minute intervals between 2007 and 2017 with the exception of 2009 and 2010 when the data were recorded at an hourly time step.

To determine the soil water storage in terms of soil water column height (SWCH in units of length) at each soil moisture sampling location (horizontally), the average soil moisture was calculated from all depths within the soil column and then multiplied by depth to bedrock at each sensor location using the depth to bedrock data from Pelletier and Rasmussen [28]. From the daily SWCH time series at each sensor location (horizontally), the sub-catchment scale daily SWCH values were estimated by computing the average of the SWCH values from all sensors within that sub-catchment. Subsequently, the MGC-scale daily SWCH values were estimated as the area-weighted average daily SWCH values of each sub-catchment (Figure 2D).

2.3. Long-Term Stable Water Isotopes in Precipitation and Streamflow

2.3.1. Precipitation Stable Water Isotope Data

Bulk precipitation samples for stable water isotope analysis were collected at multiple stations within and near MGC [15]. At the Fern Valley, Granite and Schist bulk precipitation sampling stations, two collectors were installed at each station and samples were collected every 5 to 7 days (see Heidbüchel, et al. [34] and Lyon, et al. [38] for more details). At the Mt. Lemmon and MG-Weir ISCO autosampling stations, daily bulk precipitation samples were collected. The Mt. Lemmon station was active between 2007 and 2014, but the data density decreased after 2012 for logistical reasons. The MG-Weir ISCO autosampling station has been active since 2009, but the data density similarly decreased after 2012.

In order to avoid evaporative enrichment, mineral oil was added to all water samples prior to 2015 [38]. After mid-2015, oil was not added at the Fern Valley, Schist, Granite, Mt. Lemmon and ISCO locations, in order to avoid interference with solute analyses. Samples were transported to the University of Arizona Department of Hydrology and Atmospheric Sciences laboratory within hours of collection and stored in a refrigerator at 4 °C; there is no indication of evaporative enrichment in samples collected without oil. Samples were analyzed using a DLT-100 laser spectrometer, Los Gatos Research, Inc., model #908-0008 with analytical precision (1σ) of ±0.37 ‰ and ±0.12 ‰ for δ²H and δ¹⁸O, respectively [38]. All measurements were standardized relative to international standards VSMOW and SLAP [39]. As in previous work [34], the catchment-scale time series of precipitation stable water isotopes (Figure 2A) consists of unweighted averages of observations at all five sites.

2.3.2. Stream Water Stable Water Isotope data

Stream water samples were collected using an autosampler installed at the MG-Weir site from 2006 through 2012 and stream water grab samples after 2012 (Figure 2B). Water collected in the autosampler was protected from evaporation using mineral oil. Grab samples were capped tightly on collection to prevent evaporation. The samples were analyzed using the DLT-100 laser spectrometer (see above).

3. Computational Methods

3.1. Spectral Analysis of A Time Series

3.1.1. Transformation of A Time Series

The hydrologic fluxes and storage time series data at MGC are not normally distributed (see Section S1 for more details). As a result, the data were transformed using a variance equalizing scaled arcsine hyperbolic transformation function (Equation (2)) for time series A(t) with unequal variances [19]:

$${\mathrm{A}}_{\mathrm{tra}}\left(\mathrm{t}\right)={\log }_{\mathrm{e}}\left(\frac{\mathrm{A}\left(\mathrm{t}\right)}{{\mathrm{A}}_{\mathrm{ref}}}+\sqrt{1+{\left(\frac{\mathrm{A}\left(\mathrm{t}\right)}{{\mathrm{A}}_{\mathrm{ref}}}\right)}^2}\right),$$

(2)

The variance equalizing property of the transformed time series or A_tra(t) produced by Equation (2) depends on the reference value ${\mathrm{A}}_{\mathrm{ref}}$ used in the transformation [40]. Accordingly, a sensitivity analysis of ${\mathrm{A}}_{\mathrm{ref}}$ showed that the median of a time series A(t) produced better results in terms of equalizing variance at MGC (more details in Section S2) (Figure 3).

3.1.2. The Weighted Wavelet Transform (WWT) Method for Spectral Analysis of A Time Series

The fractal scaling behavior of each time series was determined using power spectral analysis [3,19,41,42,43]. The weighted wavelet transform (WWT) method of Kirchner and Neal [19] was used due to its robustness and ability to estimate the power spectra of a time series with uneven time steps (see Section S3). Initially, the WWT method was used to estimate the period (λ)-dependent Fourier coefficients, i.e., the coefficients a (λ) and b (λ):

$$\mathrm{A}\left(\mathrm{t}\right)=\mathsf{\vartheta }+\mathrm{a}\left(\mathsf{\lambda }\right)\times \cos \left(\frac{2\mathsf{\pi }}{\mathsf{\lambda }}\mathrm{t}\right)+\mathrm{b}\left(\mathsf{\lambda }\right)\times \sin \left(\frac{2\mathsf{\pi }}{\mathsf{\lambda }}\mathrm{t}\right),$$

(3)

where $\mathsf{\vartheta }$ is a constant term for time series data. Subsequently, the power spectrum or ${\mathsf{\varphi }}_{\mathrm{A},\mathrm{wlet}}\left(\mathsf{\lambda }\right)$ of each time series was estimated by:

$${\mathsf{\varphi }}_{\mathrm{A},\mathrm{wlet}}\left(\mathsf{\lambda }\right)=0.5\left({\mathrm{a}}_{\mathrm{wlet}}^2+{\mathrm{b}}_{\mathrm{wlet}}^2\right)\times \mathrm{T}\times \frac{{\mathrm{N}}_{\mathrm{eff}}}{\mathrm{N}},$$

(4)

where T is the total time period covered by the N observations in the time series, and N_eff is the effective number of points. In Equation (4), ${\mathsf{\varphi }}_{\mathrm{A},\mathrm{wlet}}$, a_wlet and b_wlet are the estimated spectral power and Fourier coefficients for each wavelet center for the variable A(t). The period T was then divided into 100 equal intervals and values of ${\mathrm{a}}_{\mathrm{wlet}}\sqrt{\mathrm{T}\frac{{\mathrm{N}}_{\mathrm{eff}}}{2\mathrm{N}}}$, ${\mathrm{b}}_{\mathrm{wlet}}\sqrt{\mathrm{T}\frac{{\mathrm{N}}_{\mathrm{eff}}}{2\mathrm{N}}}$ and ${\mathsf{\varphi }}_{\mathrm{A},\mathrm{wlet}}$ were estimated for the center of each interval [19]. Weighted-mean values of ${\mathrm{a}}_{\mathrm{wlet}}\sqrt{\mathrm{T}\frac{{\mathrm{N}}_{\mathrm{eff}}}{2\mathrm{N}}}$, ${\mathrm{b}}_{\mathrm{wlet}}\sqrt{\mathrm{T}\frac{{\mathrm{N}}_{\mathrm{eff}}}{2\mathrm{N}}}$, and ${\mathsf{\varphi }}_{\mathrm{A},\mathrm{wlet}}$ were computed using N_eff-3, or the effective degrees of freedom, as the weighting factor. The weighted means were used to determine the final Fourier coefficients a($\mathsf{\lambda })$ and b($\mathsf{\lambda })$ and the power spectra ${\mathsf{\varphi }}_{\mathrm{A}}\left(\mathsf{\lambda }\right)$; a weight of zero was used if N_eff was less than 3 [19]. All spectral and cospectral analysis calculations were made using MATLAB® (The MathWorks, Inc., Natick, MA, USA) and an adapted version of the code of [19].

The power spectrum of a time series is commonly distorted and can be significantly inflated by spectral aliasing. Spectral aliasing describes a phenomenon by which spurious spectral power for frequencies higher than the Nyquist frequency, or the highest resolvable frequency in a time series [44], appears as spectral power associated with frequencies lower than the Nyquist frequency. For this reason, Kirchner [45] suggested using a power-law-based power spectrum filter with power-law scaling as a function of period. Accordingly, all power spectra in this work were spectrally smoothed (see Section S3) and alias-filtered following [45].

3.1.3. Statistical Hypothesis Testing

In order to evaluate the statistical significance of the relationship between a power spectrum and period for a transient variable, the correlation coefficient R² between power spectrum amplitude and period in a log-log space was tested against a correlation coefficient of zero, using the p-value as a criterion of significance at the 95% confidence level. In all cases, the null hypothesis was zero correlation and the alternative hypothesis was that the correlation coefficient was significantly different from zero.

3.2. Spectral Analysis of The Time Series of A Conservative Tracer

The spectral analysis of conservative tracer (e.g., δ¹⁸O) time series data was similar to that described in Section 3.1, but no transformation was applied to the data [20]. Spectral analysis was also avoided during time periods with significant data gaps (see Section S3) when spectral analysis would violate the Nyquist sampling theorem ([20,45]); filling data gaps by interpolation may lead to inflation in the power spectra [46]. As a result, spectral analysis of conservative tracers in precipitation and streamflow was restricted to the period between water years 2008 and 2012 due to large data gaps after water year 2012 (Figure 2A,B).

3.3. Estimation of The Period-Dependent Hydrologic Phase Index (HPI)

To determine the phase relationship of two time series A(t) and B(t), we generally followed the mathematical development of Kirchner [47] for a seasonal tracer cycle but considered multiple periods. In this method, the Fourier coefficients a_i and b_i are estimated for any two variables C_i and C_j by:

$${\mathrm{C}}_{\mathrm{i}}\left(\mathrm{t}\right)={\mathrm{a}}_{\mathrm{i}}\cos \left(\frac{2\mathsf{\pi }}{\mathsf{\lambda }}\mathrm{t}\right)+{\mathrm{b}}_{\mathrm{i}}\sin \left(\frac{2\mathsf{\pi }}{\mathsf{\lambda }}\mathrm{t}\right)+{\mathsf{\vartheta }}_{\mathrm{i}}={\mathrm{A}}_{\mathrm{i}}\sin \left(\frac{2\mathsf{\pi }}{\mathsf{\lambda }}\mathrm{t}-{\mathsf{\psi }}_{\mathrm{i}}\right)+{\mathsf{\vartheta }}_{\mathrm{i}},$$

(5)

$${\mathrm{C}}_{\mathrm{j}}\left(\mathrm{t}\right)={\mathrm{a}}_{\mathrm{j}}\cos \left(\frac{2\mathsf{\pi }}{\mathsf{\lambda }}\mathrm{t}\right)+{\mathrm{b}}_{\mathrm{j}}\sin \left(\frac{2\mathsf{\pi }}{\mathsf{\lambda }}\mathrm{t}\right)+{\mathsf{\vartheta }}_{\mathrm{j}}={\mathrm{A}}_{\mathrm{j}}\sin \left(\frac{2\mathsf{\pi }}{\mathsf{\lambda }}\mathrm{t}-{\mathsf{\psi }}_{\mathrm{j}}\right)+{\mathsf{\vartheta }}_{\mathrm{j}},$$

(6)

where ${\mathsf{\vartheta }}_{\mathrm{i}}$ and ${\mathsf{\vartheta }}_{\mathrm{j}}$ are constant terms and ${\mathsf{\psi }}_{\mathrm{i}}$ and ${\mathsf{\psi }}_{\mathrm{j}}$ are period ($\mathsf{\lambda }$)-dependent phase angles. The amplitude ratio and phase difference between variables i and j are estimated by:

$$\frac{\sqrt{{\mathrm{a}}_{\mathrm{j}}^2+{\mathrm{b}}_{\mathrm{j}}^2}}{\sqrt{{\mathrm{a}}_{\mathrm{i}}^2+{\mathrm{b}}_{\mathrm{i}}^2}}=\sqrt{\frac{{\mathsf{\varphi }}_{\mathrm{j}}^2\left(\mathsf{\lambda }\right)}{{\mathsf{\varphi }}_{\mathrm{i}}^2\left(\mathsf{\lambda }\right)}},$$

(7)

$${\mathsf{\psi }}_{\mathrm{j}}-{\mathsf{\psi }}_{\mathrm{i}}=\mathrm{phase}\left(\mathrm{i},\mathrm{j}\right)=\left(\frac{180}{\mathsf{\pi }}\right)\left\{\arctan \left(\frac{-{\mathrm{a}}_{\mathrm{j}}\left(\mathsf{\lambda }\right)}{{\mathrm{b}}_{\mathrm{j}}\left(\mathsf{\lambda }\right)}\right)-\arctan \left(\frac{-{\mathrm{a}}_{\mathrm{i}}\left(\mathsf{\lambda }\right)}{{\mathrm{b}}_{\mathrm{i}}\left(\mathsf{\lambda }\right)}\right)\right\},$$

(8)

The phase angle ($\arctan \left(\frac{-{\mathrm{a}}_{\mathrm{j}}\left(\mathsf{\lambda }\right)}{{\mathrm{b}}_{\mathrm{j}}\left(\mathsf{\lambda }\right)}\right)$ or $\arctan \left(\frac{-{\mathrm{a}}_{\mathrm{i}}\left(\mathsf{\lambda }\right)}{{\mathrm{b}}_{\mathrm{i}}\left(\mathsf{\lambda }\right)}\right)$ for variables j and i, respectively), was calculated using the MATLAB^(R) function atan2 that considers all four quadrants while estimating the inverse tangent of a function with known vertical and horizontal components, e.g., terms a and b in Equation (8) [8]. If a phase angle ${\mathsf{\psi }}_{\mathrm{j}}-{\mathsf{\psi }}_{\mathrm{i}}$ was found to be less than zero, then it was expressed as (360 + ${\mathsf{\psi }}_{\mathrm{j}}-{\mathsf{\psi }}_{\mathrm{i}})$^o. Given that our benchmarking tests suggested the potential for incorrect results when transformed data were used to determine phase differences (more details in Section S4), only un-transformed time series data, i.e., time series data without application of Equation (2), were used for this procedure.

To aid in determining if two time series are in phase or out of phase, we propose a new non-dimensional metric called Hydrologic Phase Index (HPI). Period-dependent HPI is calculated using Equation (9) below and its value ranges from –1 (when two time series are perfectly out of phase) to 1 (when two time series are perfectly in phase):

$$\mathrm{HPI}\left(\mathsf{\lambda }\right)=\cos \left(\mathrm{Phase}\left(\mathrm{i},\mathrm{j}\right)\left(\mathsf{\lambda }\right)\right),$$

(9)

where Phase (i,j) is the phase difference (see Equation (8) above) and cos is the standard cosine function. The significance of HPI values is as follows:

• HPI = 1: Two time series are perfectly in phase at some period or a they have a phase lag of 0 radians (fundamental phase lag), 2π, 4 π, … etc.;
• HPI = 0.5: Two time series are phase lagged by (π/3) radians (fundamental phase lag), 5(π/3), 7(π/3), 11(π/3), 13(π/3), … etc.;
• HPI = 0: Two time series are phase lagged by (π/2) radians (fundamental phase lag), 3(π/2), 5(π/2), 7(π/3), … etc.;
• HPI = −0.5: Two time series are phase lagged by (2π/3) radians (fundamental phase lag), 2(2π/3), 4(2π/3), 5(2π/3),7(2π/3), 8(2π/3), … etc.;
• HPI = −1: Two time series are phase lagged by π radians (fundamental phase lag), 3π, 5π, 7π, 9π, … etc.

The HPI is highly variable at any period, and it was therefore smoothed using a Gaussian smoothing window with a frequency width of 0.05 log-units, which is similar to the value used by [19] for power spectrum smoothing.

3.4. Examination of Self-Averaging Behavior in A Time Series

The average value of a time series at any given point can be calculated using moving sampling windows of width τ. Self-averaging behavior arises where the calculated averages converge on a stable value as τ is increased. To assess the potential for self-averaging behavior of the water balance components at MGC, this study applied a MATLAB® (The MathWorks, Inc., Natick, MA, USA)-adapted version of the code “RMS_averages.c” [19]. In this method, a time series is divided into τ equal windows and for each window, four different local averages are computed in a staggered pattern by shifting the window by τ/4 at each time step. Subsequently, the root mean square (RMS) differences between successive local averages are computed as a function of the timescale (τ’s). If self-averaging behavior is evident, then a slope of –0.5 is expected on a plot of linear RMS differences vs. log τ (see also Section S5), following the central limits theorem [19]. If the time series is non-stationary and non-self-averaging, then a zero or positive slope would be expected. A test for self-averaging behavior using the RMS difference method for synthetic stationary and non-stationary time series is found in Section S5 of this work (see also [19], Supporting Information).

4. Results

4.1. Power Spectra

Fractal scaling behavior was present in the power spectrum of each component of the water balance equation at high frequencies (low periods) but deteriorated at longer periods (Figure 4). In particular:

(a)

The precipitation power spectrum (Figure 4A) demonstrated statistically significant fractal scaling between spectral power and period for periods up to ~0.04 years (14 days); in this range its spectrum resembled theoretical red noise (see Supplementary Information). For periods greater than 0.04 years, the slope of the power spectrum as a function of period was approximately zero. For periods greater than 0.3 years, spectral power fluctuated greatly, and the power spectrum resembled white noise (Figure 4A).

(b)

The streamflow power spectrum showed a statistically significant relationship with period up to a period of 0.5 years. Note that the fractal slope for periods ranging from 0.16 years to 0.50 years was ~1.5 times higher than the fractal slope for the periods shorter than 0.16 years. There was no relationship for periods longer than 0.5 years (Figure 4B).

(c)

Similar to streamflow, significant fractal scaling between the AET power spectrum and period was only evident at periods less than 0.5 years (Figure 4C). For periods between 0.16 years and 0.5 years, the fractal slope was ~3.3 times higher than the fractal slope for periods shorter than 0.16 years. Although spectral analysis of the AET time series was based on a post-processed version of the observed time series data (Section 2.2.3; Knowles, et al. [35]), the power spectra for the processed AET and unprocessed AET time series for common periods were similar (see Section S6 in SI).

(d)

Fractal scaling behavior of the daily ΔS_soil was very similar to that of precipitation (Figure 4D) insofar as both showed a statistically significant relationship between spectral power and period up to a period of ~0.04 years, a white noise-like power spectrum for periods longer than ~0.04 years, and reduced spectral variability for periods between 0.04 years and 0.3 years relative to longer periods.

4.2. Spectral Analysis of δ¹⁸O in Precipitation and Streamflow

There was statistically significant (p << 0.001) evidence of fractal scaling in the precipitation-δ¹⁸O power spectrum for periods less than 0.2 years, but only white noise thereafter (Figure 5). In contrast, the stream water-δ¹⁸O power spectrum showed statistically significant (p << 0.001) fractal scaling at all periods, which was not the case in the in the stream water flux at periods longer than 0.5 years (Figure 4B). Across all periods, the fractal slope of 1.07 for δ¹⁸O in stream water was within a previously reported range (1.05 ± 0.15) for various conservative and non-conservative solutes [19].

4.3. Phase Relationships as A Function of Component Period

Phase shifts between precipitation and the other water balance components evolved with component period and did not assume a constant value (Figure 6). At sub-seasonal scales (periods less than 0.1 years), the HPI spectrum suggested a variety of phase relationships and significant impact from spectral smoothing (see Supplementary Information). Consequently, we restricted the current analysis to consider only phase relationships at periods longer than 0.1 years. The HPI curves at 5 to 10 year periods (Figure 6) do not yield meaningful information at MGC because the period length approaches the length of the entire dataset.

The HPI values between precipitation and streamflow were positive at all periods (Figure 6A) and the minimum HPI (0.07) occurred at a period of 0.7 years. Seasonally, i.e., for period = 0.5 year (because there are two wet seasons per year), streamflow lagged precipitation by 63° (~1.1 radians). Between precipitation and AET, HPI was slightly negative at periods between 0.1 years and 0.6 years, except for periods between 0.14 and ~0.2 years when HPI was positive. At periods greater than 0.6 years, HPI generally increased with period. Seasonally, AET lagged precipitation by 86° (~1.5 radians). HPI between precipitation and daily ΔS_soil was generally positive. Seasonally, ΔS_soil lagged precipitation by 53° (~0.9 radians).

4.4. Self-Averaging Time Scales

For all water balance fluxes, log-log plots of RMS difference as a function of averaging time interval τ showed a slope of –0.5 (the expected slope for a self-averaging time series) only past a threshold τ value, which was different for every component (Figure 7). The specific threshold τ values were ~0.04 years for precipitation, ~0.2 years for streamflow, and 0.35 years for AET. For daily ΔS_soil, the log-log slope of RMS difference versus averaging interval τ was –0.5 for τ values up to 0.2 years and was steeper thereafter. Interestingly, the daily soil water storage (S_soil, as distinct from the water balance component daily change in storage, ΔS_soil, see Equation (1)) showed no evidence of self-averaging behavior at any τ value (see Supplementary Information). For S_soil and for larger values of averaging time scales (up to 2.8 years), the slope was near zero.

5. Discussion

Power spectral analysis with respect to period indicated two principal behaviors in the daily water-balance components at MGC (Figure 4; Figure 8). Red noise, or significant fractal slopes at shorter periods, was observed at periods less than 0.04 years for precipitation (P) and ΔS_soil, and periods less than 0.50 years for stream runoff and AET. At these periods, the spectrum for each variable was positively correlated with period. In contrast, white noise, or random variations in the amplitude of fluctuations, was observed for P and ΔS_soil components at periods longer than 0.04 years and for Q and AET components at periods longer than 0.5 years. A detailed characterization of period-dependent spectral behavior and the corresponding phase relationships between specific water balance components is provided below.

5.1. Temporal Variability of Various Water Balance Components and Their Phase Relationships

The fractal part of the precipitation spectrum was associated with wetter and drier intervals of about 7 days each (7 days is half the 0.04-year period length), which may represent a characteristic length of such intervals at MGC. Up to a period of 14 days (0.04 years), periodic fluctuations cumulatively sampled shorter time intervals, which resulted in more precipitation in longer periods, and the subsequent red noise behavior. At periods longer than 14 days, periodic fluctuations generally sampled only one 7-day wet interval; hence there was no cumulative precipitation increase and a zero fractal slope. We observed high-amplitude spectral power at the seasonal period of ~0.5 years that corresponded to the cumulative precipitation of the whole monsoon or winter season. Patterns of soil wetting closely followed precipitation events with damping and phase shifts at shorter periods resulting from advection-dispersion effects in soils that are not typically saturated (Figure 4). At longer periods (more than 0.3 years), the ΔS_soil power spectra showed less damping and was less out of phase with precipitation, resulting from saturation of the soil. Phase shifts and damping of ΔS_soil relative to precipitation reflect the tendency of soils to retain water even during dry periods.

Both AET and runoff (Q) are typically generated from soil water at MGC [15,16]. As a result, additional damping and/or phase-shifts in Q and AET relative to precipitation or ΔS_soil at short periods (less than 0.5 years) was due to retention of soil water during unsaturated conditions. Lower phase shifts for Q relative to AET corresponded to direct runoff of precipitation; there was little phase-shift of Q relative to P at short periods. We therefore ascribe transmission of high-amplitude white noise from the P signal to rapid flow of water through saturated soils (AET and Q) or runoff (Q). In contrast to the ΔS_soil spectrum, the S_soil spectrum showed a strongly persistent behavior, i.e., significant correlation between power spectrum and period with fractal slope >1 at lower periods (up to 0.5 years), and a weakly persistent behavior that was not statistically significant at higher periods (see Supplemental Information). In this way, soil water storage at MGC functioned as a low-pass filter at periods up to 0.3 years with much less filtering of the input signal at periods greater than 0.5 years. Filtering behavior of the soil reservoir at MGC was thus a direct result of soil water content.

5.2. Observed Power Spectra within the Context of Previous Work

The observed fractal scaling behavior of precipitation, AET, streamflow, and the daily change in soil water storage (Figure 4) was broadly similar to the results of previous work that included only one or two of water balance components in their analysis. For example, Kirchner, et al. [7], Weedon, et al. [8], and Milly and Wetherald [9] reported that the streamflow power spectra were damped relative to the precipitation spectra up to threshold periods, beyond which the two power spectra were more closely aligned. Our results for precipitation, streamflow, and AET are also similar to those in large river basins located in the Russian Federation, where there are strong fractal relationships at shorter periods, but white-noise spectra at longer periods [10]. The transition from fractal scaling to white noise behavior for the daily change in soil water storage at MGC occurred at a period of ~0.04 years, compared with estimates of 1 to 2 years for mid-latitude soils, calculated from a model-derived time series [11]. Insofar as the time series of all water balance components at MGC were based on long-term observations and considered together, this work advances and extends understanding of observed hydrological fractal scaling behavior for a mountain headwater environment.

5.3. Dynamic Phase Relationships between the Various Water Balance Components

Our results support the existence of dynamic phase relationships between the water balance components at MGC (Section 4.3). For example, Q and P were close to being perfectly in phase at 2–3 years, but out of phase by 90° (¼ period, or 3 months) at ~1 year (Figure 6), and Q lagged P by 60° at the time scale of rain events or wet intervals (7 to 60 days). Likewise, P and AET were in phase at a ~3 year period, but ~90 out of phase at periods less than 0.7 year, and P and ΔS_soil, were out of phase by 60° (⅙ cycle, or ~ 18 days) and 45° (⅛ cycle, or ~ 3 months) for periods of 0.3 years and 2 years, respectively. Accordingly, we caution against fixed or period-invariant interpretation of phase relationships that could obscure the dynamic nature of catchment-scale hydrologic fluxes and water partitioning.

The dimensionless HPI, proposed in Section 4.3, provides a convenient way of visualizing dynamic phase shifts, and should be applicable to hydrologic model calibration and other disciplines in which time-series data are converted to power spectra. Additionally, this metric could be useful in the development, calibration, and diagnostics of hydrologic models. For example, a recent study involving cross-spectral analysis between observed and simulated discharge time series and assessment of performance of several hydrologic models in the Thames basin, UK [8] used a phase range of −180° to 180°, both extremes corresponding to an out-of-phase relationship. In the same situation, if the proposed HPI metric is used, then a range of HPI from −1 to +1 for perfectly out-of-phase and perfectly in-phase relationships, respectively, is arguably a clearer way of expressing phase difference.

5.4. Self-Averaging Behavior

The threshold τ values in Figure 7 correspond to the inception of self-averaging behavior, and not to values of τ for which an adequate estimate of each average could be obtained. For example, the observed seasonal and annual variability of P are such that ten or more years of data would be required to obtain an accurate estimate of the average precipitation at MGC. The specific threshold τ value for P was ~0.04 years, which corresponded to the period at which there was a slope change in the P power spectrum (Figure 4A), and self-averaging behavior for τ > 0.04 years was consistent with this behavior (Section 5.1). Approximately similar fractal slope breaks and threshold τ behavior was observed for Q (Figure 4B and Figure 7B), AET (Figure 4C and Figure 7C) and ΔS_soil (Figure 4D and Figure 7D). Thus, in a broader sense, self-averaging of the various water balance components was related to their respective power spectral characteristics.

5.5. Marshall Gulch Catchment (MGC) as A Fractal Filter

MGC behaved as a fractal filter insofar as it converted random white noise in the P input signal into a 1/f stream water output “signal” [19]. This result was confirmed by the δ¹⁸O time series data where δ¹⁸O in precipitation showed a white noise power spectrum at periods longer than 0.2 years, but the δ¹⁸O power spectrum in stream water was significantly correlated with period at all periods (Section 4.2; Figure 5). Further, the fractal slope for the stream water δ¹⁸O was close to 1, which is indicative of 1/f scaling behavior [7,18,19,20,45]. Note that 1/f filtering behavior in stream water chemistry is shown to stem from subsurface advection and dispersion (both in space and time) of spatially-distributed random chemical inputs to a catchment [19]. Previous work at MGC determined that the δ¹⁸O in catchment-scale precipitation was dominated by temporal rather than spatial variations, and that both spatial and temporal dispersion occurred in the subsurface where flow-path lengths and transit times differed across a range of scales [48]. As a result, observed fractal filtering behavior at MGC likely results from a combination of these processes.

The 1/f filtering behavior for a conservative tracer at MGC was similar to the observed behavior for conservative tracers (e.g., chloride) in other humid catchments and sites located in coastal settings in North America and Europe [7,18,19,20]. This behavior may reflect a dominant contribution of soil water to streamflow at these sites, which has been shown at MGC for both streamflow and seasonal vegetation water [15,16], and is supported by the short residence time of soil water storage at MGC [48] and in humid and coastal settings [20]. Overall, subsurface water storage was the dominant control on the fractal filtering behavior of hydrological fluxes and tracer concentration in outflows (e.g., streamflow and evapotranspiration) relative to the power sprectra of precipitation input. The current study therefore reinforces the viewpoint that a thorough understanding of subsurface water storage is necessary to understand hydrologic functioning across climatic, geologic and geomorphic settings [29]. Similarly, groundwater stored in fractured and/or weathered bedrock is also likely to factor into fractal filtering of catchments precipitation inputs ([7] and references therein).

6. Conclusions

We determined power spectra from catchment-scale time-series data of precipitation (P), changes in soil water storage (ΔS_soil), actual evapotranspiration (AET) and streamflow (Q) using the weighted wavelet transform method. The P power spectra showed fractal scaling for periods up to 0.04 years and white noise thereafter. Otherwise, the power spectra of ΔS_soil, AET and Q were damped at short periods and phase-shifted relative to the P input signal. In contrast to previous findings, phase-shifts among the water-balance components were dynamic functions of period. White noise in the P-δ¹⁸O input signal was filtered by the catchment and emerged as a 1/f flicker signal in the Q-δ¹⁸O output signal. Similar fractal filtering was observed between the P and Q time series data where white noise at periods of 0.04 to 0.3 years in the P power spectrum emerged as a fractal signal in the Q power spectrum. Catchment-scale fractal filtering resulted from both spatially and temporally varying tracer and flux inputs and subsurface advection and dispersion effects arising from different flow paths and wetting properties of soil. All water balance components showed self-averaging behavior, but at time scales that differed among the various components. Thus, this study constrains the interactions of the various components of the hydrologic cycle in a mountainous catchment; such interactions are fundamental to further study of the hydrologic and hydro-bio-geochemical functioning of such catchments.