Journal of Hydrology xxx (2014) xxx–xxx Contents lists available at ScienceDirect Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol Spatial and temporal structure within moisture measurements of a stormwater control system Ruben Kertesz ⇑, Lee Rhea, Daniel J. Murray Jr. U.S. EPA, National Risk Management Research Laboratory, 26 West Martin Luther King Drive, Cincinnati, OH 45268, United States a r t i c l e i n f o Article history: Available online xxxx Keywords: Pervious concrete Soil Sensor Moisture Stormwater Infiltration s u m m a r y This study develops novel geostatistical methods to investigate the spatial relationship between individual soil moisture sensors placed within native soil and #57 crushed stone aggregate subbase. The subbase sensors are beneath a 0.06 ha (0.15 acre) pervious concrete parking lot in Cincinnati, OH, USA. The parking lot treats runon from a 0.198 ha (0.49 acre) asphalt area. A geostatistical characterization of moisture (measured as permittivity) in the subbase beneath pervious concrete indicates that significant spatial correlation is either not present or only present at very short distances (<2.5 m). A twostage para-statistical model relating antecedent storm moisture to apparent pervious concrete infiltration was developed to identify temporal trends in the data and to detect the clogging processes with relatively simple parameterization. The results suggest that either the placement of the sensors is not sufficient to detect clogging or that clogging is not problematic for the study period. Suggestions are provided to improve future research installations, based upon the findings here. Subbase moisture analysis results are compared with native soil moisture results. Seasonal trends are more pronounced in the native soil than in the subbase. The statistical analyses are applicable to multiple Storm Control Measures (SCM), Best Management Practices (BMP), agriculture, and soil environments. Other studies can determine the statistical power of their sensor installation using the methods applied here, which are flexible enough for multiple applications. Furthermore, data reduction methods presented serve to easily elucidate short-term moisture responses due to rainfall. A quantile response pattern is provided for sensors installed in both subbase and soil. Published by Elsevier B.V. 1. Introduction Moisture measurement in soils has been an important agricultural pursuit since the early 20th century (Ayers et al., 1943) and is commonly used in viticulture (Matese et al., 2009; Lunt et al., 2005). Various sensors have been developed to measure moisture and sensor accuracy has been consistently improved for better performance in soils (Seyfried and Murdock, 2001; Kelleners et al., 2005), but their use in very well drained granular media is a relatively novel application of the technology and few studies have examined their performance in such environments (e.g. Spittlehouse, 2000; Arguedas et al., 2007). This is a nascent issue in the stormwater field because most previous green infrastructure studies, where moisture measurements have been made, tended to focus on bioretention (Houdeshel and Pomeroy, 2010; Brown and Hunt, 2011; Darner and Dumouchelle, 2011; Aravena and Dussaillant, 2009). However, recent applications include use in ⇑ Corresponding author. Tel.: +1 513 569 7428. E-mail address: kertesz.ruben@epa.gov (R. Kertesz). green roof media (Palla et al., 2011; Sun et al., 2013) and beneath pervious concrete (Brown and Borst, 2013; Stander et al. 2013). Previous research in stone aggregate or coarse media has been performed using largely different sampling network design. Brown and Borst (2013) observed the down gradient migration of clogging along a roadway curb of interlocking concrete pavers at a site in Louisville, KY. The site had a very narrow and long flowpath which was investigated using a geometric transect in the well-defined flowpath. Research by Palla et al. (2011) and Sun et al. (2013) focused on the drying and wetting patterns of green roofs with comparatively thin media depths. Stander et al. (2013) described the arrangement of water content reflectometers beneath permeable pavement. The network design involved sensor placed near the pavement edge for different permeable pavement types. All of the aforementioned involve a unique spatial arrangement of sensors, but to the knowledge of the authors little if any information is available regarding optimization of sensor placement with respect to statistical concerns such as sample independence or spatial correlation. Studies focused on overall effects would logically optimize sensor placement by using http://dx.doi.org/10.1016/j.jhydrol.2014.02.011 0022-1694/Published by Elsevier B.V. Please cite this article in press as: Kertesz, R., et al. Spatial and temporal structure within moisture measurements of a stormwater control system. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2014.02.011 2 R. Kertesz et al. / Journal of Hydrology xxx (2014) xxx–xxx adequate sensor spacing to ensure spatial independence, and thereby maximize statistical power. Conversely, studies focused on aspects such as media uniformity might optimally space sensors by ensuring they were within the distance over which spatial correlation was expected to extend. Dobriyal et al. (2012) have documented multiple indirect moisture measurement technologies that suffer from limited spatial representation in soil. On the scale of most field investigations, moisture sensors are often thought of as point samples. The present study presents information pertinent to the determination of the representativeness of individual sensor measurements; if spatial correlation within a system is present, it could be reasonable to treat the readings from individual sensors as representative of some volume of soil within the range of spatial correlation, rather than merely a number of point-wise observations. Among the differing technologies available to measure subsurface moisture, the technology applied in this study is transmission line oscillation (TLO) water content reflectometery [Campbell Science CS-650]. TLO technology is based upon the well-described research by Topp et al. (1980), who formulated an equation to calculate volumetric water content using permittivity. Advancements have been made in moisture sensing technology since then (Campbell and Anderson, 1998), and Kelleners et al. (2005) are among the latest to describe TLO moisture measurement methods. Many other sensor technologies have been developed to measure real-time moisture as well, such as: tensiometry (Coolong et al., 2012); time domain reflectometry (TDR) (Malicki et al., 1996; Robinson et al., 2003); time domain transmission (Blonquist et al., 2005); and frequency domain reflectometry (Bandaranayake et al., 2007). TLO water content reflectometers were selected for the present study because they are manufactured for the datalogging platform deployed at the site and were readily available. The authors do not intend to form comparisons between technologies or make any claim as to the (dis)advantages of a particular moisture sensing technology. Rather, the methods presented here are intended to be cross-applicable. The promise of gathering data from sensors in a stone aggregate subbase beneath pervious concrete (PC), is in being able to use the timing of response to determine patterns in the treatment system. This may be relatively easier to identify in bioretention systems, which have a distinct wetting front that may take between minutes and days to migrate through the soil (Hsieh and Davis, 2005; Dussaillant et al., 2005; Shuster et al., 2007; Stander et al., 2010; Brown and Hunt, 2011), subject to differences in the gradation of the bedding layer and subsurface media (Robinson and Friedman, 2001). However, the movement of water through the open-graded subbase (stone aggregate) beneath pervious concrete is almost instantaneous (Chopra et al., 2010; Stander et al., 2013; Brown and Borst, 2013). Selective placement of sensors is critical to answering many research questions. Examples include: determining the spatial progression of clogging; comparing rainfallresponse timing through permeable pavement to the response in soil; identifying drainage areas that contribute a disproportionate volume of runoff to one area of the SCM; or identifying localized flooding in a SCM. The present study evaluates the utility of moisture sensors to ascertain spatial and temporal moisture trends in a stone aggregate subbase beneath a functioning pervious concrete (PC) parking lot under varying weather conditions. It also evaluates wet weather responses in native soil. The objectives of this study are to: develop methods to determine the presence of spatial correlation in sensor response (as permittivity) monitored in stone aggregate beneath pervious concrete; determine the utility of these data for monitoring pavement clogging using the established spatial network of sensors; contrast readings with those made in native soil; and identify the ramifications for instrumentation and monitoring frequency requirements in future studies. 2. Materials and methods 2.1. Site description The site is located at latitude 39°80 22.9400 N, longitude 84°250 46.700 W, in the northwest region of the Cincinnati Public School’s Clark Montessori High School. The subject PC parking lot is part of a stormwater treatment train to reduce discharge to combined sewers. A 0.198 ha (0.49 acre) asphalt drainage area contributes flow to both the subject 0.061 ha (0.15 acre) PC, and a 0.028 ha (0.07 acre) bioretention area (BA). The dimensions of stormwater control system are shown in Fig. 1. Flowpath lines were generated using a ground-based LIDAR survey and are documented in Fig. 1 as well. Overflow drains through a raised outlet at the far west end of the BA if ponding occurs in the BA. The site is part of a larger combined effort by the Metropolitan Sewer District of Greater Cincinnati (MSD) to evaluate green infrastructure (GI) as a tool to reduce combined sewer overflows. 2.2. Instrumentation and data management Many sensor technologies on the market calculate dielectric (electrical) permittivity as a surrogate for volumetric water content. The Campbell Scientific sensors deployed here calculate dielectric permittivity as a function of Period and electrical conductivity. Stander et al. (2013) recommend the use of permittivity instead of volumetric water content when installing moisture sensors in stone aggregate because of difficulty calibrating the sensors in this type of material (Stenger et al., 2005). The permittivity data recorded in the current investigation are affected by occasional erratic responses in electrical conductivity readings, perhaps due to interference from electrical currents in subsurface conduits. Baselining serves to mitigate the effect of that behavior. Statistical methods shown in this manuscript are presented for permittivity but are also applicable to soil moisture. Sensors are deployed along three east–west transects parallel to the longitudinal axis of the PC, designated Upstream (U), Downstream (D), and Native (N) (Fig. 1). The U- and D-transects are located 152 mm (6 in.) and 457 mm (18 in.) from the upgradient pavement edge, respectively. The N-transect is located along the centerline of the associated drainage swale. The sensors within the U- and D-transects are emplaced in a well-drained #57-stone subbase beneath the pervious concrete (152 mm above the native soil). The sensors within the N-transect are emplaced in densified soil (denoted native soil) approximately 150 mm into the native soil. The native soil type is a stiff to very stiff mixture of brown, lean clay, with traces of organics, rock, and brick fragments. Table 2 shows physical properties of the pervious concrete. The sensors are geometrically spaced along the longitudinal axis of transects. The moisture sensors were installed during the construction process. The U- and D-transect sensors were installed by first placing them and ancillary communications cabling on the soil surface. Next, a protective cylinder was placed around them before lifting sensors, adding media to the inside of the cylinder to a depth of 152 mm (6 in.), and laying sensors on the media bed. The sensors were then covered with media until flush with the media depth surrounding the cylinder. Lastly, the cylinder was lifted up and out, leaving the sensors installed in the matrix before pervious concrete was poured (Fig. 2). The N-transect sensors were installed by pushing the sensor rods (tines) horizontally into the native soil, perpendicular to the BA, after the BA was excavated. Please cite this article in press as: Kertesz, R., et al. Spatial and temporal structure within moisture measurements of a stormwater control system. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2014.02.011 R. Kertesz et al. / Journal of Hydrology xxx (2014) xxx–xxx 3 Fig. 1. Horizontal (A) and vertical (B) distribution of sensors at the site. Spaces are numbered from 1 to 22 (East to West) along the 60 m (198 ft.) transect. Numerical notes for part B are as follows. (1) 305 mm (12 in.) open graded #57 stone; (2) filter fabric lapped onto gravel surface; (3) flush concrete curb; (4) subgrade 92% proctor; (5) 152 mm (6 in.) pervious concrete. Flowlines are superimposed upon (A), colored by drainage location. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 1 Various parameters that were logged at each sensor. Parameter Unit Volumetric water content Electric conductivity Temperature Permittivity Period Voltage ratio m3/m3 dS/m °C Unitless lS V/V Table 2 Physical properties of pervious concrete. Parameter Density (kg/m3) Average (n = 8) Standard deviation 2186 82 Mix design test Value (various units) Design strength (psi) Water/cement ratio Cementitious (% volume) Aggregate (% volume) Design air (% volume) Slump Aggregate type 3000 0.31 11 58 20 0 #8 Gravel All cabling was placed in conduit before installation and connected to a centralized CR-1000 datalogger (Campbell Scientific, Inc., Logan, Utah). Six parameters were logged via serial digital interface protocol (Table 1) and telemetrically uploaded to a centralized server. Rainfall data measured after October 2011 were recorded at 1 min intervals using a 152 mm (6 in.) diameter tipping bucket rain gauge with a 0.1 mm gradation (R.M Young Company, Traverse City, Michigan). Five minute rainfall data were provided by MSD’s Rainvieux system for dates before this time, a service of Vieux Inc. (Norman, OK). The service technology is discussed by Gourley and Vieux (2005). In order to maintain a Fig. 2. Illustration of sensor installation in porous media before filling cylinder and surrounding infiltration gallery. uniform dataset, the relatively brief record (8/23/11–10/21/11) of 5-min Rainvieux data were linearly apportioned to 1-min increments to match the substantial record of 1-min data collected for the remainder of the investigation. Data were uploaded from dataloggers to a remote server using a wide area network connection and then loaded into a database with an appropriate format for data reduction and statistical analysis using Base SASÒ v.9.3 (SAS Institute Inc., Cary, NC, USA). Data reduction and statistical analyses were performed using the 1-min timestep data recorded at the site. The statistical methods are timestep agnostic. 2.3. Data reduction Data reduction for the present study was performed in three stages. In the first stage, the response data were baselined for each individual sensor to remove low-frequency effects. Such effects would obfuscate high-frequency, storm-related responses. In the Please cite this article in press as: Kertesz, R., et al. Spatial and temporal structure within moisture measurements of a stormwater control system. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2014.02.011 4 R. Kertesz et al. / Journal of Hydrology xxx (2014) xxx–xxx second stage, transfer functions were developed for individual sensors that described the time-convolved relationship between precipitation and sensor response. In the third stage, the ensemble of the adjusted parameters for all sensors was analyzed in a single model to account for spatiotemporal correlations in the data. The model was later used to discern whether some effect only ascribable to sensor distance from the upgradient edge of the pavement (i.e. clogging) was present. This step-wise procedure was intended to minimize the influences of sensor-specific variations on the overall analysis for the presence of clogging. 2.3.1. Baselining Baselining is an important step in rainfall response analysis. Sensor baseline readings vary considerably, obfuscating stormspecific responses which are often lagged or convolved in time. Baselining data allows the researcher to remove low frequency (or long-term, e.g. weeks to months) components in the signal to elucidate any (relatively) simple high-frequency (short term) relationships between individual rainfall events and sensor response. Initially, raw permittivity was plotted for each sensor against timestamps and manually reviewed for anomalies such as inexplicable jumps in baseline (dry weather) values, apparent drifts in sensor sensitivity (due to factors such as instrument aging), and confounding effects related to weather (such as air temperature and potential evaporation). Pronounced variation in the instrument baselines was observed, so all time series were re-baselined using: quadratic functions of the temperature of the monitored media; date stamps; and a linear function of a sine-cosine pair with a period of one year. The error function minimized was the sum of squared estimation for positive errors (overestimates) and 1/100 this value for negative errors (underestimates). The ratio was applied because underestimates were not as pronounced as overestimates. Fitting was done using SAS proc model. The calculated baseline was subtracted from the raw instrument responses to avail comparison of data between instruments, with respect to time. 2.3.2. Rainfall-response deconvolution The baselined rainfall-permittivity response can be deconvolved into components using a transfer function. Differentials in component timestamps between sensors then provide insight into changes in sensor response with respect to time as well as space, such as clogging. Different transfer function methods were evaluated to model the relationship between precipitation and infiltration through pervious concrete into the subbase. There are parameterizations of such processes that impose restrictions on the form of the transfer function based on numerous assumptions, but these methods generally suffer from an inability to accurately reflect unknown systematic deviations from their underlying assumptions, such as unanticipated bimodal or trimodal responses. Another class of transfer function estimation relies on using the linear covariances or related quantities between the response of interest and lags of precipitation (such as least squares, linear programming, matrix decompositions, and others). While these methods generally have no problems with ‘‘accurately’’ capturing atypical behaviors in transfer functions, they do have problems with stability and over-dispersion. Using individually calculated covariances between the response and various lags of precipitation will generally provide a plausible transfer function. However, it will often be over-dispersed due to autocorrelation across lags, and therefore indicate responses start before the precipitation fell. Techniques based on solving a linear system using statistically-based methods do correct for over-dispersion but suffer from stability problems because of the difficulties with inverting a design matrix with column-wise collinearity. These problems are generally unavoidable if lags of a response variable are used because the response is usually highly autocorrelated. The problem is typically less severe if lags of precipitation are used, but still remains if the observation interval is shorter than the duration of an individual storm. These problems were circumvented herein by using a novel approach to estimate the transfer functions needed for the model using the primary right eigenvector (v) from singular value decompositions of the matrix direct-products of the baselined moisture readings (Q), with a Toeplitz matrix comprised of columns of the precipitation observations offset by increasing lags (Ptoep): Q  Ptoep, where  denotes a matrix direct-product. Eigenvectors were similarly calculated, including those for Date ðd  PtoepÞ, and a sinusoid with a wavelength of one year (sin(y) Ptoep þ cosðyÞ  PtopÞ. The date eigenvector is intended to capture any temporal change in responses attributable to clogging. The sinusoidal eigenvector is intended to capture seasonal interactions such as temperature-dependence. Initial results showed that the calculated primary eigenvectors had similar forms to published rainfall-runoff transfer functions that were developed using more typical techniques such as ordinary least squares. In comparison to the aforementioned transfer functions, however, the righthand-side eigenvectors appeared to be insignificantly affected by either of the over-dispersion or stability issues. This is possibly due to the matrix rotations inherent in generating the eigenvectors. The primary right hand side eigenvectors were therefore adopted for the present study as the transfer functions relating rainfall and sensor response. Precipitation was transformed using the eigenvectors to remove the temporal offsets between precipitation and instrument responses. 2.3.3. Geospatial modeling The primary objective of the geostatistical analysis is to discern whether the data from each sensor can be treated as spatially independent and therefore amenable to classical statistical analyses, or if, in addition to temporal autocorrelation, spatial correlation exists that must be accounted for in statistical tests. Overlay plots of precipitation and moisture indicate that the near-saturation population is associated with rainfall, so a rainfall intensity of 0.1 mm per minute was used to separate the dry-versus wet-weather data for the geostatistical analysis. Normal probability plots (Fig. 3) were used to further divide moisture readings from U- and D-subbase sensors into four populations. Based upon the plot shown in Fig. 3, the sensor responses in the subbase appear to be from the following four distributions: (1) a horizontal limb corresponding to instrument baseline or very dry weather (permittivity less than 1); (2) a gradually rising limb corresponding to dry weather with residual subsurface moisture (permittivity between 1 and 4); (3) a vertical limb corresponding to a set of readings that bridge the residual subsurface moisture and subsurface saturation (permittivity between 4 and 30); and (4) a relatively flat limb corresponding to subsurface measurements approaching saturation (permittivity >30). Note that some events do not reach permittivity measurements greater than 5 at either upstream or downstream locations (e.g. 4/21/2012). The probability plots indicate that the first two populations include over 99% of the data. In contrast, a very small proportion of the data on the probability plot falls on the nearly vertical line that is the third population and few data are saturated (fourth population). A fraction of the third population could be a transitional group where water surrounding the stone aggregate drains rapidly by gravity. However, the authors find that it accurately represents the immediate short-term response to a storm event. It is the third, near-vertical, segment between residual moisture and quasi-saturation that is interesting in the context of clogging analysis. Hence, analyses were based on the data from the population associated with immediate response to precipitation. Please cite this article in press as: Kertesz, R., et al. Spatial and temporal structure within moisture measurements of a stormwater control system. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2014.02.011 5 R. Kertesz et al. / Journal of Hydrology xxx (2014) xxx–xxx Fig. 3. Normal quantile plots for base-lined moisture (as permittivity) at typical stone aggregate (Upgradient = U, Downgradient = D) and soil (Natural Soil = N) transects. Sensor responses belonging to a normally distributed group should fall on a straight line. Maximum saturation is typically reached at a permittivity of 40. date effects such as periodicity or offsets in response times between sensors. All of these models have limitations in their ability to fit semivariograms. They are thus suboptimal for the purpose of determining whether there is any identifiable semivariance or lack of independence in the measured data (Isaaks and Srivastava, 1989). Consequently, the authors sought a generalized semivariogram model more flexible than any of the aforementioned models; the model would thereby provide greater statistical power to ascertain whether there is lack of independence between moisture readings from different sensors. We adopted the cumulative density function of the continuous Weibull distribution as the semivariogram model. This function includes a shape parameter as a continuous variable that can be tuned to not only mimic the aforementioned variogram models, but nearly any intermediate between them. Expressing this function in terms of semivariance as a function of lag: t The geostatistical analysis is limited by the spatial distribution of the data. The moisture data for soil are clearly from a different population than the moisture data for stone aggregate, so an initial division of the data was made by media type. Consequently the analysis is limited in spatial dimensions, as the sensors used in the study were organized into three parallel transects with the two in the subbase being separated by no more than one foot, and each having all sensors placed at the same elevation above the native soil. These transects are spaced in this fashion to capture the progress of an expected ‘‘clogging front’’, anticipated to move from the upgradient pavement edge toward the center of the parking lot. The semivariogram calculations for the moisture in the subbase therefore include: (1) the longitudinal separation between sensors along the (parallel) axes of the upgradient and downgradient transects; and (2) the transverse separation between the upgradient and downgradient transects. Calculations do not include a depth component, but they do allow for anisotropy between the longitudinal and transverse axes relative to the parallel transects. This determination was made after assessing whether a semivariance model could be fit to the semivariogram calculated from the moisture data, the results of which indicated a significant relationship between semivariance and lag. The semivariogram was calculated as one-half the average squared difference between contemporaneous moisture readings at each possible pairing of sensors. The spatial distance between each pairing of sensors was used as the corresponding lag. Typical semivariogram models include components for short-range variance unattributable to lag (‘‘nugget’’), and a trend of increasing semivariance with increasing lag. Common basic models include linear, spherical, exponential, and Gaussian. Many models are curvilinear to accommodate a semivariogram that approaches a maximum semivariance (the ‘‘sill’’) asymptotically; this is equivalent to the ‘‘true’’ variance of the data in the absence of cross-correlation between sensors. The linear model, however, optimally fits only a constant trend in semivariance with lag (and optionally a nugget). A fit to the semivariogram by this model generally indicates that either there is a trend in the data that should be corrected for prior to semivariance calculations, or that the maximum separation of the data points was insufficient to fully characterize the spatial structure of the data. In contrast, the spherical model optimally fits only semivariograms that reach independence within a finite lag. It often overestimates the rate of increase in semivariance at relatively short lags. The exponential model, and its cousin the gaussian model, differ from the linear and spherical models in that they both approach the sill asymptotically; the Gaussian model includes an inflection that allows for persistent semivariance at short lags. None of these basic models accommo- DðxÞ ¼ k  keð1½xu=sÞ þ nugget ð1Þ where x (m) is lag; u (m) is x-offset; s = scale parameter; t = shape parameter (recall (x  u)/s is the standard normal deviate for Gaussian data), and k = scale factor normalizing D(x) to interval [0, 1]. This is a modification of the more familiar 3-parameter expression for this function: YðxÞ ¼ 1  eð½xu=sÞ t ð2Þ Our modifications provide a monotonic function that can include a nugget and the maximum value of the response variable (D, or Y in the equations above). The three parameter Weibull was chosen as a basis because it allows an offset on the x-axis to control for variation in sensor response due to longitudinal position. Differential clogging between the upgradient and downgradient transects could be manifested as very high anisotropy. This model was therefore fit to the semivariogram data by iteratively recalculating the lags as Euclidean distances with a varying constant multiplied by the (transverse) distance component. 2.3.4. Clogging After first baselining and deconcolving the data, then developing a model to account for (auto)correlation with respect to space and time, data were analyzed to investigate the occurrence of clogging. Time-series plots of moisture sensor permittivity readings within the subbase indicated that the number of apparent infiltration events, or instances of maximum instrument response, were less frequent during the latter part of the study relative to the earlier portion of the study. Some instruments might have been changing in sensitivity as the study progressed. However, precipitation also appeared to accumulate at a slower rate as the study progressed (Fig. 4), raising the possibility that the decrease in the frequency of infiltration events could have been merely due to decreased rainfall. Therefore, a statistical model was employed to differentiate such nuisance effects on subbase moisture readings. The clogging model was designed to discern whether sensor responses were differentiable by some factor that might indicate clogging rather than by measurable climatic conditions or sensor aging (drift). We assumed that sensor group was a proxy for proximity to the expected clogging front; the front was assumed to progress from the upgradient pavement edge toward the center of the pavement. Sensor response and precipitation were related in the model through the previously calculated transfer functions for each sensor (time-invariant temporal convolutions). Nuisance effects included in the model were: a smoothly changing instrument response through time (a quadratic ‘‘main effect’’ of time and time-squared); the (now convolved) prior 30 min of precipitation, not changing as a function of time (a ‘‘main effect’’ Please cite this article in press as: Kertesz, R., et al. Spatial and temporal structure within moisture measurements of a stormwater control system. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2014.02.011 6 R. Kertesz et al. / Journal of Hydrology xxx (2014) xxx–xxx Fig. 4. Precipitation intensity and accumulation throughout the study. A linear regression of cumulated precipitation with date accentuates the relatively steady decrease in measured precipitation as the study progresses. DT is the date and timestamp of measured precipitation. of antecedent precipitation); and the prior 30 min of precipitation, changing as a linear function of time (an ‘‘interaction effect’’ between precipitation and time). Quadratic functions captured any change in instrument sensitivity through time (instrument drift). Linear effects were used to model time-invariant and timevariant responses to convolved precipitation. Our initial concern was that quadratic and linear models are potentially too constrained to model such effects without biases if they exhibit behaviors such as logarithmic decays or stepresponses. We found however that, when used in conjunction with a first-order autoregressive (AR1) term, they were sufficient to contain both nuisance effects and autocorrelation. This finding is supported by the non-significance (P < 0.1) of the slope of the linear regression of residuals on their first lag. We therefore assumed that the residuals were not problematically autocorrelated, indicating that standard errors of the model fixed effects should not be biased and thereby render invalid statistical tests. The significance of each term in the model was evaluated using F-tests for overall fixed effects and random effects, and t-tests for differentiating individual levels within fixed effects. 3. Results and discussion 3.1. Spatiotemporal description of data Subsurface permittivity responses differed markedly between those placed in the stone aggregate subbase and those placed in soil. In contrast with sensors placed in aggregate, sensors in soil did not exhibit sharp increases in permittivity immediately after or during storm events (short term), as evidenced by the differences shown in Fig. 3 as well as between Figs. S1 and S4. Sensors in soil also exhibited greater seasonality (e.g. Fig. 5, upgradient stone aggregate sensor V2CU versus soil sensor V2CN). While multiple years of data were not available, seasonality is evidenced by changes in event frequency that are congruent with seasonal periods (Fig. 4). Detail plots illustrating short-term sensor responses are not included in the main text for brevity and because it is extremely difficult to visually separate additive effects comprising signals. Supplemental Fig. S1, however, shows upstream and downstream pairs for selected storm events at sites V2, V5, and V10. Both upstream and downstream V2 sensors reach saturation more often than V5, which, itself, reaches saturation more often than V10. During smaller events (i.e. 4/21/2012), V5 has the highest response. At V2 and V5, the downstream response is muted, compared to upstream. The opposite is true at V10. Site V10 also shows a more linear rise to peak and decay compared to other locations. Fig. 5. Examples of long-term sub-pavement aggregate subbase (A) and sub-trench soil (B) moisture readings as permittivity. Rainfall (mm) is also plotted on the same scale. The same linearity is observed for volumetric water content (not shown). This indicates that there may be a unique drainage phenomenon at those downstream locations. However, as shown in Fig. S2 (storm responses for additional upstream locations), V22CU has a similarly linear signal. Like Fig. S1, Fig. S2 shows decreasing saturation frequency with increasing location number. The sensor responses at V1CU are saturated during more events than at location V18CU and there are no instances of saturation at site V22 for the time period shown. The apparent lack of consistent difference in response between subbase sites may be simply due to the difficulty of visually decomposing signals into their constituent parts, or due to a number of factors that are difficult to control such as non-uniformity of the bedding material, sensor divergence through long-term drift, non-uniform clogging, construction irregularities, etc. In fact, observations made during construction included a distinct difference in gravel bed depth between locations V1 and locations V5 through V22. There is a 180 mm difference in depth to surface between V1 and V5 (Fig. S3). This would certainly lead to higher magnitude peaks and a sharper rise to peak at shallow beds because the bed would get saturated more easily than locations with a deeper bed depth under storms of the same intensity. Such localized phenomena may diminish the strength of geospatial trends between sensors. 3.2. Spatial data structure Semivariograms were constructed from upstream and downstream location data using the responses for events of intensities greater than 0.1 mm min1. Semivariograms initially constructed using raw data indicated an unbounded trend in increasing semivariance with increasing lag; when using baselined data and a planar site-wide trend surface, the trend disappeared. Iterative fitting of our semivariogram model indicated that, if a zero x-intercept is assumed, the range of spatial correlation would have to be very short to accommodate the observed data (Fig. 6). Fit characteristics are tabulated in Table 3. The range is visually estimated to be between 2.5 m and 6 m but there is no significant support for the presence of any spatial correlation structure in the data. If the Please cite this article in press as: Kertesz, R., et al. Spatial and temporal structure within moisture measurements of a stormwater control system. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2014.02.011 7 R. Kertesz et al. / Journal of Hydrology xxx (2014) xxx–xxx generalized semivariogram model is not constrained (not shown), then it approximates a flat line with no significance (model fit P > 0.1). Note that many of the high semivariance values are associated with V2 Upstream (V2CU). Even when V2CU is excluded, an estimated lag of 6 m would be needed for data to be independent (when forced through x = 0). As indicated during the discussion of spatiotemporal trends, the statistical analysis shows that the sensors, as installed, do not provide enough information to show a spatial correlation structure. There is not enough information to determine if general pavement response could be determined using a couple of sensors or if there would need to be multiple sensors, separated by a minimum distance (e.g. 6 m). The pavement performance monitoring method is dictated by whether or not the researcher thinks that there should be theoretically zero variation between collocated sensors or not. 3.2.1. Clogging Differences in relatively short-term response characteristics were apparent between subbase sensors; responses were differentiated by sensor group, which could be attributed to differential clogging along the expected migration of a clogging front (Table 4). The magnitude of the group-wise difference is small, in part, due to the varying scales of the constituent variables in the model, which were not standardized. The practice of non-standardization or normalization is in conformance with the opinions of some statisticians and in nonconformance with the opinions of others. Regardless of its magnitude, the detection of group-wise difference supports the concept of using the presented method or variations of it to detect clogging. While the vast majority of the variation in the data is attributable to the (random) first-order autoregressive terms associated with lagged permittivity (Estimate > 0.97), significant amounts of information are also captured by the precipitation effect (Estimate 0.027). Effects of drift or seasonality through time (Time) are not significant with respect to individual sensors. Seasonality is likely absorbed by the autoregressive (or moving average) parameter in the model. 3.2.2. Comparison of native soil response to subbase response While it is of little use to apply the previously discussed statistical clogging analysis method to the sensors deployed in native soil, the process of baselining data revealed interesting qualitative differences between permittivity measurements in native soil and stone aggregate. Fig. 1 illustrates the locations of native soil Fig. 6. Variograms for permittivity in stone aggregate subbase during precipitation intensity at or above 0.1 mm per minute. The variogram function is the 3-parameter cumulative density function (CDF) for the Weibull distribution. A straight- line-fit (green, dashed) is provided for visual reference. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 3 Fit parameters for Weibull-CDF based variogram model and comparison to a linear model for spatial correlation of permittivity in crushed stone subbase beneath porous pavement. Neither model is statistically significant (P = 0.48, 0.43). Models were statistically differentiable from each other (P = 0.02). Parameter Value Weibull CDF fit x-Offset Nugget (y-Offset) Range (correlation length) Sill (true variance) Shape Scale Anisotropy coefficient 0.0 0.75 20928.0 143.0 0.18 15137.0 2.66 E12 Value Statistic Model fit evaluations SSE R-sq R-sq adjusted P-value Statistic Model comparison F-stat P-value Weibull-CDF Linear 49006 0.017 0.015 0.48 Value 48834 0.02 0.01 0.43 5.56 0.02 sensors along the BA’s east-west transect. Baselines differ between native soil and aggregate, indicating increased soil moisture retention in the soil (Fig. 5). This is supported by geotechnical analysis of soil cores which found that the first 2.3 m of soil consists of brown lean clay, trace organics, and rock and brick fragments. The soil, in contrast to upstream and downstream sites, has a muted shortterm response (Fig. 3) and much of the change in moisture occurs over an extended period of time. As described previously, seasonal trends are more pronounced in the native soil than in the aggregate beneath the PC (Fig. 5). It is reasonable to expect a stronger influence of seasonality on sensors in native soil because temperature has a stronger influence on two factors, permittivity and electrical conductivity, when there is more water present in the surrounding media. As temperature decreases, permittivity increases but bulk electrical conductivity decreases. The temperature dependence of both permittivity and conductivity are large. In comparing the water content of clay and gravel after gravitational water drains, there is much more water held by clay than gravel. Results also show that the short term response to precipitation varies by sensor location. Supplemental Fig. S7 shows that the magnitude of baselined response to rainfall broadly increases from V2N to V22N. While the trend is affected by aberrations in permittivity, particularly at V22CN, a brief analysis of Period data (which were largely unaffected by electrical conductivity) supports the findings and indicates that V22CN maintains a high saturation level. However, both Period and permittivity data show that site V5CN maintains a permanently high signal as well. The sensor produced readings during storm events that were higher than even a saturated clay soil (Bittelli et al., 2008). Perhaps there is persistent ponding at site V5CN, or perhaps localized soil chemistry has affected the sensor. Barring site V5CN, the lag effect is congruent with the physical construction of the facility. Due to the slope of the bioswale towards the west and due to the purposeful blockage of the underdrain outflow, upgradient free water is expected to collect in the underdrain and exfiltrate back out of the underdrain downgradient. One major difference between the responses in the native soil and those in the subbase are that the native soil sensors exhibit unusual perturbations that are difficult to explain, while upstream and downstream subbase locations along the same vector do not. Please cite this article in press as: Kertesz, R., et al. Spatial and temporal structure within moisture measurements of a stormwater control system. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2014.02.011 8 R. Kertesz et al. / Journal of Hydrology xxx (2014) xxx–xxx Table 4 Statistical summary for model of baselined moisture sensor readings (permittivity). P-values less than or equal to 0.01 are considered significant. ‘Group’ is a proxy for proximity to an expected clogging front. Sensor response and precipitation were related through temporal convolution prior to running the statistical model. Effect type Effect-level-specific effects Group Fixed effects Group Group Sensor (Group) Sensor (Group) Sensor (Group) Sensor (Group) Sensor (Group) Sensor (Group) Sensor (Group) Sensor (Group) Sensor (Group) Time (Date2) Time (Date) Temperature (Celcius2) Temperature (Celsius) Precipitation Intercept Random effects Lagged e Lagged e Lagged e Lagged e Lagged e Lagged e Lagged e Lagged e Lagged e D U D D D U U U U U U Sensor V2CD V5CD V10CD V1CU V2CU V5CU V10CU V18CU V22CU V2CD V5CD V10CD V1CU V2CU V5CU V10CU V18CU V22CU Overall effects Estimate Std. Err. t-Value P-value F-value ProbF 0.0013 0 0.0028 0.0054 0 0.0076 0.0004 0.0053 0 0.0005 0.0013 9.30E19 3.10E09 0 0.0003 0.0267 2.5877 0.0003 na 0.0003 0.0003 na 0.0003 0.0003 0.0003 na 0.0003 0.0003 na na na 0 0.0009 1.0872 4.13 na 9.91 18.14 na 26.75 1.60 18.43 na 1.72 4.18 na na na 13.99 30.75 2.38 <.01 Na <.01 <.01 Na <.01 0.11 <.01 Na 0.09 <.01 Na Na Na <.01 <.01 0.02 19.65 <.01 280.15 <.01 10.12 <.01 243.46 <.01 945.78 5.67 <.01 0.02 0.9753 0.9736 0.9893 0.9724 0.9842 0.9727 0.9852 0.9847 0.9890 0.0002 0.0003 0.0004 0.0002 0.0001 0.0003 0.0003 0.0002 0.0003 4821.45 3402.63 2780.93 4985.77 8059.12 3793.06 2884.19 4603.28 2918.41 <.01 <.01 <.01 <.01 <.01 <.01 <.01 <.01 <.01 nc nc Some effects of very low magnitude are significant because they were not standardized (e.g. Time/Date). High sensor-specific autocorrelation was contained by specifying the first time-lag of individual sensor responses as the only random effect. na: parameter estimate not available (may be zero-level automatically selected by the software). nc: not calculated but apparently significant. e: dimensionless permittivity. These perturbations deviate sharply from the strong seasonal pattern described previously. Supplemental Figs. S5–S7 show longterm data series for a number of upstream, downstream, and native sites. Site V5CN shows a distinct rise and then sudden drop in permittivity response in June, 2012. V22CN shows a continuously elevated permittivity response throughout the first part of 2012. It is unknown whether this was a measurement of real hydrologic phenomenon or not. It is also important to note that the baselined response is designed to normalize results for better comparison. However, there is a possibility that it obfuscates real phenomena (e.g. V22CN). For example, February 2013 readings of many native soil sites show a sharp drop in permittivity that is not observed in baselined data. This is a common effect of freezing water but native soil temperatures at the sample sites were measured to drop to 2–3 °C, just above freezing. 3.2.3. Contribution to future research The implications of these findings on future research are that, while it appears that there is no clear geospatial trend in the sensors placed beneath the pavement at the test site, the methods developed here should be applied to other investigations. The Weibull function is flexible, such that it can be used to analyze other multi-sensor datasets for geospatial trends without modification. Particular emphasis should be placed on spacing sensors along more densely populated transects to determine if it results in better correlation and lower variance between sensors when installed in coarse media. Improvements to the current study include: placing sensors in a gravel bed of succifient depth; designing the pavement surface with zero transverse slope; placing sensors along a geometrically increasing distance from the pavement edge using multiple lag distances; and installing sensors in duplicate. Lastly, the results suggests that maximum power is achieved by separating treatments from each other by >6 m. Another very interesting tool that can be applied to other research methods includes the observations made in Fig. 3, where the normal quantile plots indicate distinct phases of moisture response. Researchers can produce such plots to quickly expose values of interest and to focus on specific components of the drying/wetting cycle. Similarly, the method used for baselining data may be valuable for application to other datasets to account for long-term sensor drift and seasonal effects and is a useful first step in data reduction. The methods used to detect clogging are a particularly novel approach to the problem and could be applied to other installations, again with a different spatial arrangement of sensors. In particular, the presented statistical methods are unique in that they avoid issues of overdispersion while also remaining stable. They allow modeling of both time-variant and time-invariant responses to precipitation and can be adjusted to include coupled parameters not investigated in this project. All of the methods above can be broadly applied to analyze other types of moisture sensors and sensor response variables as well as other applications. Clogging or sealing investigations may be an interesting component of floodplain analysis in order to determine if effective infiltration is changing over time. 4. Conclusions Methods are developed to extract useful spatial and temporal information from moisture sensors embedded in soil and crushed stone aggregate. A model removes seasonal trends and offsets. A simple and intelligent data visualization and reduction technique is developed. A geostatistical model is used to investigate sensor (in)dependence. The methods are appropriate for application to various moisture sensor technologies and soil conditions. Please cite this article in press as: Kertesz, R., et al. Spatial and temporal structure within moisture measurements of a stormwater control system. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2014.02.011 R. Kertesz et al. / Journal of Hydrology xxx (2014) xxx–xxx It is apparent that, due to the lack of typical soil properties (capillary action, matric suction, etc), moisture sensors respond differently in crushed stone aggregate (in subsurface stormwater retention/detention areas) than in soil. An investigation of the spatial distribution of sensors installed beneath PC suggests that sensor independence is achieved within a short lag distance (6 m). There appears to be more than 1 phase of response in most sensors buried beneath PC. Statistical analysis suggests that there is not enough information to identify an upstream-downstream clogging pattern during the investigation period. A better arrangement of sensors should lead to a more conclusive outcome. Disclaimer The US Environmental Protection Agency, through its Office of Research and Development, funded and managed, or partially funded and collaborated in, the research described herein. It has been subjected to the Agency’s administrative review and has been approved for external publication. Any opinions expressed in this paper are those of the author (s) and do not necessarily reflect the views of the Agency, therefore, no official endorsement should be inferred. Any mention of trade names or commercial products does not constitute endorsement or recommendation for use. Acknowledgments This research was conducted with the support of an appointment to the Research Participation Program at the National Risk Management Research Laboratory administered by the Oak Ridge Institute for Science and Education. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jhydrol.2014.02. 011. References Aravena, J., Dussaillant, A., 2009. Storm-water infiltration and focused recharge modeling with finite-volume two-dimensional richards equation: application to an experimental rain garden. J. Hydraul. Eng. 135 (12), 1073–1080. Arguedas, F.R., Lea-Cox, J.D., Ristvey, A.G., 2007. Characterizing air and water content of soilless substrates to optimize root growth. Comb. Proc. Int. Pl. Prop. Soc. 57, 701–708. Ayers, A.D., Wadleigh, C.H., Magistad, O.C., 1943. The interrelationships of salt concentration and soil moisture content with the growth of beans. Agron. J. 35, 796. Bandaranayake, W.M., Parsons, L.R., Borhan, M.S., Holeton, J.D., 2007. Performance of a capacitance-type soil water probe in a well-drained sandy soil. Soil Sci. Soc. Am. J. 71, 993–1002. Bittelli, M., Salvatorelli, F., Pisa, P.R., 2008. Correction of TDR-based soil water content measurements in conductive soils. Geoderma 143, 133–142. Blonquist Jr., J.M., Jones, S.B., Robinson, D.A., 2005. A time domain transmission sensor with TDR performance characteristics. J. Hydrol. 314, 235–245. Brown, R., Borst, M., 2013. Assessment of clogging dynamics in permeable pavement systems with time domain reflectometers. J. Environ. Eng. 139 (10), 125–1265. 9 Brown, R.A., Hunt, W.F., 2011. Impacts of media depth on effluent water quality and hydrologic performance of under-sized bioretention cells. J. Irrig. Drain. E. ASCE 137 (3), 132–143. Campbell, G.S., Anderson, R.Y., 1998. Evaluation of simple transmission line oscillators for soil moisture measurement. Comput. Electron. Agr. 20, 31–44. Chopra, M., Kakuturu, S., Ballock, C., Spence, J., Wanielista, M., 2010. Effect of rejuvenation methods on the infiltration rates of pervious concrete pavements. J. Hydrol. Eng. 15 (6), 426–433. Coolong, T., Snyder, J., Warner, R., Strang, J., Surendran, S., 2012. The relationship between soil water potential, environmental factors, and plant moisture status for poblano pepper grown using tensiometer-scheduled irrigation. Int. J. Veg. Sci. 18, 137–152. Darner, R.A., Dumouchelle, D.H., 2011. Hydraulic characteristics of low-impact development practices in northeastern Ohio, 2008–2010. U.S. Geological Survey Scientific Investigations Report 2011-5165. USGS, Virginia. Dobriyal, P., Qureshi, A., Badola, R., Hussain, S.A., 2012. A review of the methods available for estimating soil moisture and its implications for water resource management. J. Hydrol. 458–459, 110–117. Dussaillant, A.R., Cuevas, A., Potter, K.W., 2005. Rain gardens for stormwater infiltration and focused groundwater recharge: simulations for different climates. Water Sci. Technol. Water Supply 5 (3), 173–179. Gourley, J., Vieux, B.E., 2005. A method for evaluating the accuracy of quantitative precipitation estimates from a hydrologic modeling perspective. J. Hydrometerol. 6, 115–133. Houdeshel, C., Pomeroy, C., 2010. Predicting and Verifying Physiological Performance of Plants in a Bioretention System in an Arid Climate. World Environ. Water Resour. Congress 2010, 4019–4024. Hsieh, C., Davis, A., 2005. Evaluation and optimization of bioretention media for treatment of urban storm water runoff. J. Environ. Eng. 131 (11), 1521–1531. Isaaks, E.H., Srivastava, R.M., 1989. Applied Geostatistics. Oxford University Press, New York. Kelleners, T.J., Seyfried, M.S., Blonquist, J.M., Bilskie, J., Chandler, D.G., 2005. Improved interpretation of water content reflectometer measurements in soils. Soil Sci. Soc. Am. J. 69, 1684. Lunt, I.A., Hubbard, S.S., Rubin, Y., 2005. Soil moisture content estimation using ground-penetrating radar reflection data. J. Hydrol. 307, 254–269. Malicki, M.A., Plagge, R., Roth, C.H., 1996. Improving the calibration of dielectric TDR soil moisture determination taking into account the solid soil. Eur. J. Soil Sci. 47, 357–366. Matese, A., Di Gennaro, S.F., Zaldei, A., Genesio, L., Vaccari, F.P., 2009. A wireless sensor network for precision viticulture: the NAV system. Comput. Electron. Agr. 69, 51–58. Palla, A., Sansalone, J.J., Gnecco, I., Lanza, L.G., 2011. Storm water infiltration in a monitored green roof for hydrologic restoration. Water Sci. Technol. 64, 766. Robinson, D.A., Friedman, S.P., 2001. Effect of particle size distribution on the effective dielectric permittivity of saturated granular media. Water Resour. Res. 37 (1), 33–40. Robinson, D.A., Jones, S.B., Wraith, J.M., Or, D., Friedman, S.P., 2003. A review of advances in dielectric and electrical conductivity measurement in soils using time domain reflectometry. Vadose Zone J. 2, 444–475. Seyfried, M.S., Murdock, M.D., 2001. Response of a new soil water sensor to variable soil, water content, and temperature. Soil Sci. Soc. Am. J. 65, 28. Shuster, W.D., Gehring, R., Gerken, J., 2007. Prospects for enhanced groundwater recharge via infiltration of urban storm water runoff: a case study. J. Soil Water Conserv. 62 (3), 129–137. Spittlehouse, D.L., 2000. Using time domain reflectometry in stony forest soil. Can. J. Soil Sci. 80 (1), 3–11. Stander, E., Borst, M., O’Connor, T., Rowe, A., 2010. The effects of rain garden size on hydrologic performance. World Environ. Water Resour. Congr. 2010, 3018– 3027. Stander, E., Rowe, A., Borst, M., O’Connor, T., 2013. Novel use of time domain reflectometry in infiltration-based low impact development practices. J. Irrig. Drain Eng. 139 (8), 625–634. Stenger, R., Barkle, G., Burgess, C., 2005. Laboratory calibrations of water content reflectometers and their in-situ verification. Aust. J. Soil Res. 43, 607. Sun, T., Bou-Zeid, E., Wang, Z.-H., Zerba, E., Ni, G.-H., 2013. Hydrometeorological determinants of green roof performance via a vertically-resolved model for heat and water transport. Build. Environ. 60, 211–224. Topp, G.C., Davis, J.L., Annan, A.P., 1980. Electromagnetic determination of soil water content: measurements in coaxial transmission lines. Water Resour. Res. 16 (3), 574–582. Please cite this article in press as: Kertesz, R., et al. Spatial and temporal structure within moisture measurements of a stormwater control system. J. Hydrol. (2014), http://dx.doi.org/10.1016/j.jhydrol.2014.02.011