Introduction

Groundwater levels are fundamental to the practice of hydrogeology. They are necessary for flow system analysis and understanding from determining hydraulic gradients to estimating material properties to calibrating regional groundwater flow models.

From a hydrological perspective, groundwater levels and their seasonal fluctuation affect runoff generation in southern Ontario. Longer-term (5-10yr) fluctuations could also have impacts on land use development; think of constructing and maintaining buried infrastructure (e.g., sewers, water mains), or applying certain Low Impact Design (LID) strategies, etc.

For a number of good reasons, groundwater resource management frequently tends to boil down to a set of steady-state (static) targets. This immediately implies that the steady-state condition represents some long-term average state, without determining, say: How much should I expect the water table to vary from its current position?

The ORMGP

The Oak Ridges Moraine Groundwater Program reflects the intent of fifteen government agencies (municipalities and conservation authorities) to jointly manage, synthesize and disseminate water resources data, analysis, knowledge, geologic interpretations, and numerical models over a 30,000 km² area situated in south-central Ontario.

Interpretive ORMGP maps showing estimates of the depth to the water table and deeper potentiometric surfaces are gaining increased use. These maps are constructed from thousands of data points collected over many years and at various times of the year and are considered to represent an “average” state. To effectively use these maps, users must be aware of variability in the groundwater system and consider i) longer-term deviation from the measured and average conditions, ii) short-term seasonal patterns, and iii) water table sensitivity to storm and/or snow melt events.

From the program’s database, 689 good-quality, active and inactive, monitoring wells exist within the study area, collecting regular (e.g., daily, hourly or less) datalogger groundwater level data over decades. At some locations, when evaluated alongside climate data (e.g., precipitation and snow melt estimates) the water level response to events or stresses can be analysed. An example location, Alton PS 3-1, is provided below.

A 4.2m deep well located 120m away from Shaw’s Creek, Alton, ON.


Longer-term levels

Long-term seasonal variability in shallow groundwater levels has been quantified using a 12-point Generalized Additive Model (GAM) that provides confidence intervals to a fitted seasonal trend.

From 273 shallow wells (depths < 20m) analysed, the south-central Ontario water table fluctuates ±0.7 (s.d. 0.6) m relative to some “average” state. Deeper (>20m depth; n=416) wells show greater fluctuation of groundwater levels at ±1.4 (s.d. 1.6) m.

Methodology

Pattern + Noise

The intention here is to determine the expected range of variability one should expect when interpreting the ORMGP interpolated water table/potentiometric surface maps. The water table map should be thought of as a long-term average with uncertainty associated with when the measurements were made. It is a surface interpolated from 100,000s of measurements made over the past century, and constrained by mapped watercourses, water bodies, subsurface sediment permeability, etc.

Realistically, the depth to water table and potentiometric surfaces can vary naturally by a metre or even more, especially where water levels are affected by local groundwater pumping centres. Long-term trends in water levels and aquifer storage can further impact expected water levels, as this changes the stage from which the expected variability is to occur.

For now, lets skip the event-based groundwater response and break the expected variability into 2 components:

\[\text{inter-annual variability} + \text{seasonal variability}\]

Further more, lets assert:

  • seasonal variability is assumed to have an annual cyclic nature; and
  • inter-annual variability (or long-term trend) can be used to de-trend the time-series.

Generalized Additive Model

The Generalized Additive Model (GAM) is a flexible means of fitting a model to time series data. It is a generalization of the linear models we have all fit to a scatter plot while retaining the ability to efficiently compute confidence/prediction intervals using an intuitive control (Wood, 2017).

GAMs are especially useful as hydrographs contain a good deal of noisy auto-correlated data (that is, a sequence of measurements made in a short time frame tend to be correlated). For instance, it’s apparent that there is a seasonal pattern to most hydrological data in southern Ontario. Here we design the GAM with 12 smoothing spline “knots”, one for every calendar month. These monthly knots are further specified as being a cyclic regression spline by assuring that the surface remains continuous to the second derivative at a year’s end (Wood, 2017). The GAM will fit a smoothing spline that captures the underlying seasonal pattern.

The same goes for the long-term trend, which is determined using knots specified to every year. The GAM is then specified as:

\[ \texttt{h}_i = f_{yr}(\texttt{year}_i) + f_\tilde{m}(\texttt{day.of.year}_i) + e_i, \quad e_i=\phi e_{i-1}+\epsilon_i \]

where \(h_i\) is the expected/modelled water level, the “smooth function” \(f_{yr}\) is the long-term trend and \(f_\tilde{m}\) is the seasonal variability of the GAM, and \(e_i\) is the AR(1) autocorrelation function. The GAM is initially set to having 365 + \(n\) degrees of freedom, \(n\) is the number of years.

What’s important is that the above equation is linear, and so the GAM can be decomposed into long-term and seasonal patterns fitted to the data:

An example GAM analysis outputted for NVCA - Midhurst BH2 (PGMN W0000244-2).


The above four plots illustrate an example of the output for a PGMN groundwater level monitoring location (NVCA - Midhurst BH2). The bottom 2 plots illustrate the first 2 terms \(\left(f_{yr}(\texttt{year}_i) \text{ and } f_\tilde{m}(\texttt{day.of.year}_i)\right)\) of the above equation, with 90% prediction intervals shaded. The top 2 plots show the GAM vs. observed both over the period of record (left) and after being de-trended (right).

The y-axis label of the bottom 2 plots report numbers called the “effective degrees of freedom,” which relates to the minimum number of constraints needed to reproduce this GAM, simplifying the regression from the initial constraints given (i.e., \(\text{e.d.f}=3.4+9.5 \ll 365+n\).)

Processing steps

The ORMGP database currently hosts 689 locations with >34 water level logger measurements having a well screen depth and at least 4 years of data. In addition, for every location time-series:

  1. outliers were automatically removed using the ±1.5 IQR (inter-quartile range) method.
  2. water level data were re-sampled to daily averages.
  3. values are factored according to their day-of-year and year of the reported groundwater level.
  4. data are de-trended by subtracting \(f_{yr}\) from the data.
  5. 90% prediction intervals surrounding the GAM are added to the smooth function \(f_\tilde{m}\).
  6. annual maximum and minimum extent of the prediction intervals defines the anticipated range in groundwater levels one should expect at a given location in the near future. (This is given in the upper-right plot in the above figure.)

Maps

Water Level Range

hover over any circle to reveal station properties, click for more details. Full-screen available in the top-left corner

Water Table

GAM modelled seasonal groundwater level variability for shallow (<20m) wells.


Deeper system groundwater levels

GAM modelled seasonal groundwater level variability for deeper (>20m) wells.


Distributions


By formation


Summary table, groundwater range (±m) by formation
FORMATION n mean geomean median
Late Stage Glaciolacustrine-Glaciofluvial 16 0.58 0.43 0.43
Halton Till 34 0.64 0.49 0.51
Mackinaw/Oak Ridges 169 0.67 0.43 0.44
Channel - Silt 15 0.89 0.59 0.43
Channel - Sand 19 1.01 0.53 0.48
Upper Newmarket 30 0.64 0.52 0.63
Inter Newmarket Sediment 34 0.68 0.47 0.51
Lower Newmarket 66 0.72 0.55 0.65
Thorncliffe 131 1.87 1.04 1.12
Sunnybrook 21 1.09 0.63 0.49
Scarborough 80 1.92 1.11 1.08
Bedrock - Undifferentiated 72 0.88 0.64 0.64


Long-term Trend (last 10 years)

Groundwater level monitoring over the past decade has been put to the Mann-Kendall test for trend. The maps below identify wells that have exhibited a statistically significant trend (increasing or decreasing or none) change over the past 10 years.

Mann-Kendall test for trend.

  • Blue Triangles: increasing trend (p<0.05)
  • Orange Triangles: decreasing trend (p<0.05)
  • Circles: no trend

hovering on any triangle to reveal station properties, click for more details. Full-screen available in the top-left corner

Shallow (<20m) wells


Deep (>20m) wells


Distributions



References

Wood, S.N., 2017. Generalized Additive Models: An introduction with R, 2nd ed. CRC Press. 476pp.



Extras

Complete GAM formulation

\[ g(\mu_i)=\mathbf{A}_i\boldsymbol\theta + \sum_{n=1}^{p}b_{n}(t)\beta_{n} + \sum_{m=1}^{12}\tilde{b}_m(t)\beta_m + \mathbf{Z}_i\mathbf{b} \]

where \(\mathbf{A}_i\) is the \(i^\text{th}\) row of a model matrix \(\mathbf{A}\), \(\boldsymbol\theta\) is the parameter vector. \(\boldsymbol\theta\) can be leveraged to handle pumping influences; \(\tilde{b}_m(t)\) is the cyclic basis function for month \(m\) at time \(t\) multiplied by some parameter \(\beta_m\); and \(\mathbf{Z}_i\mathbf{b}\quad(\mathbf{b}\sim N(\mathbf{0},\boldsymbol\psi))\) is the error term.

Example in ggplot

cdf <- gamm(data=exampldf, Val~s(doy,bs="cc",k=12)+s(as.numeric(yr),bs="cr"), correlation=corAR1(form=~1|yr))
exampldf %>% ggplot(aes(date,Val)) +
  theme_bw() +
  geom_line(aes(colour='observed')) +
  geom_line(data=data.frame(Val = cdf$gam$fitted.values, date = exampldf$date), aes(colour='modelled')) +
  labs(x=NULL, y="water level (masl)")

print(summary(cdf$gam))
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Val ~ s(doy, bs = "cc", k = 12) + s(as.numeric(yr), bs = "cr")
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 230.41501    0.05817    3961   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                     edf Ref.df     F p-value    
## s(doy)            9.446  10.00 101.6  <2e-16 ***
## s(as.numeric(yr)) 3.330   3.33  19.6  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.791   
##   Scale est. = 0.062726  n = 5501
layout(matrix(1:2, nrow = 1))
plot(cdf$gam, scale=0)

gam.check(cdf$gam)

## 
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                      k'   edf k-index p-value    
## s(doy)            10.00  9.45    1.03    0.97    
## s(as.numeric(yr))  9.00  3.33    0.33  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
intervals(cdf$lme,which = "var-cov")
## Approximate 95% confidence intervals
## 
##  Random Effects:
##   Level: g 
##                   lower         est.       upper
## sd(Xr - 1) 0.0003575794 0.0008996196 0.002287185
##   Level: g.0 
##                     lower        est.      upper
## sd(Xr.0 - 1) 0.0001811636 0.001539592 0.01378532
## 
##  Correlation structure:
##         lower      est.     upper
## Phi 0.9976847 0.9988042 0.9993826
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Within-group standard error:
##     lower      est.     upper 
## 0.1800430 0.2504520 0.3483956
cdf$gam$sig2/cdf$gam$sp
##            s(doy) s(as.numeric(yr)) 
##      0.0008996196      0.0015395917

simultaneous interval for a penalized spline in a fitted Generalized additive model (GAM)



see also

Groundwater and Flooding: investigating the causal nature of groundwater flooding through simulation.