In jkennel/waterlevel: Well Water Level Analysis Tools

knitr::opts_chunk$set(
  collapse = TRUE,
  fig.ext = 'png',
  comment = "#>"
)

Introduction

This vignette covers the static barometric/loading efficiency (BE/LE) calculations in the waterlevels package. We will use a dataset collected from a fractured sandstone site (called transducer) as an example which has 1.5 month of data collected every 2 minutes and columns called datetime, wl, et, and baro which corresponds to the date and time of the analysis, the water level pressure from a non-vented submerged transducer, the synthetic Earth tide, and the atmospheric pressure from a non-vented transducer. We use the raw pressure measurements and the raw barometric pressure data (i.e. there is no reason to convert to an estimated water level for these analyses).

Static barometric efficiency calculations are simple to apply but it is important to consider their applicability on a well by well basis.

In particular, these methods assume:

• That there is no time-lag in the response between the water levels and the barometric pressure. This response is often called the fully confined response. If the monitoring frequency is large in comparison to the time lag this assumption may be relaxed somewhat.

We start by loading the packages and applicable data for the loading efficiency calculation.

library(waterlevel)
library(earthtide)
library(data.table)
library(plotly)
library(recipes)
library(rlang)
library(purrr)
library(dplyr)

data(transducer)

First step: Plot the data

par(mar = c(2, 5, 0.5, 0.5))
plot(baro~datetime, transducer[seq(1, nrow(transducer), 10)], type = 'l', 
     las = 1, lwd = 2, col = '#BC7C56',
     xlab = '', ylab = 'baro (dbar)')

plot(wl~datetime, transducer[seq(1, nrow(transducer), 10)], type = 'l', 
     las = 1, lwd = 2, col = '#5696BC',
     xlab = '', ylab = 'wl (dbar)')

Remove a trend

It is often important to remove trends from the data prior to analysis. We will remove the trend in the water level data and the barometric pressure data using linear regression. Some methods are more robust for handling data with complicated trends, but for this dataset, and these analyses a linear trend should suffice. We can then replot the values. Now that we have removed the trend we can plot the water level pressure and the atmospheric pressure on the same figure for coparison.

data(transducer)
transducer[, wl_detrend := lm(wl~datetime)$residuals]
transducer[, baro_detrend := lm(baro~datetime)$residuals]

par(mar = c(2, 5, 0.5, 0.5))
plot(baro_detrend~datetime, transducer[seq(1, nrow(transducer), 10)], 
     type = 'l', las = 1,
     col = '#BC7C56',  lwd = 2,
     xlab = '', ylab = 'baro detrend (dbar)')

points(wl_detrend~datetime, transducer[seq(1, nrow(transducer), 10)], 
       type = 'l', las = 1, 
       col = '#5696BC', lwd = 2,
       xlab = '', ylab = 'wl detrend (dbar)')

The two responses look very similar. At this scale it appears that the water level response is a smoothed version of the barometric response with a slightly diminished magnitude. Let's explore this more fully.

Calculate Barometric/Loading Efficiency

There are many ways to calculate the barometric/loading efficiency. I have tried to include as many methods as possible in the waterlevel package and if you know of a method I missed let me know and I will do my best to add it.

Least squares

This method is just a simple wrapper for least squares regression added to be consistent with the other functions in the package. The inverse column is either TRUE or FALSE and deserves a few words here. The value should be set to TRUE to signify that if the independent variable (atmospheric pressure) increases the dependent variable decreases (dep). This is typically the case for vented transducers and water level obtained via a water level tape. The result of this calculation will be a barometric efficiency. For non-vented transducers when atmospheric pressure increases the submerge pressure reading also increases. The result of this calculation is the loading efficiency.

Summary:

• Loading efficiency, when inverse == FALSE (for non-vented transducers)

• Barometric efficiency, when inverse == TRUE (for vented transducers, manual measurements)

be_ls <- be_least_squares(dat = transducer, 
                          dep = 'wl_detrend', 
                          ind = 'baro_detrend',
                          inverse = FALSE,
                          return_model = FALSE)

# This value is the BE/LE
be_ls

par(mar = c(5, 5, 0.5, 0.5))
plot(wl_detrend~baro_detrend, transducer[seq(1, nrow(transducer), 10)],
     pch = 20, las = 1,
     cex = 0.05, col = '#000000',
     asp = 1)
grid()

If you are interested in the residuals from the model fit or the regression results you may specify return_model = TRUE. In this method the loading efficiency is r round(be_ls, 2).

Least squares using first differences

In the presence of background trends running the regression on first differences can minimize their effect on the regression results. To apply this method we have one more parameter called lag_space which determines how many samples are between difference calculations. This value is important when monitoring frequency is high. Commonly BE/LE calculates the differences based on data collected at 1 hour increments. So for this example where we collected 2 minute data that translates to a lag_space of 30

be_lsd <- be_least_squares_diff(dat = transducer, 
                                dep = 'wl_detrend', 
                                ind = 'baro_detrend',
                                lag_space = 30,
                                inverse = FALSE,
                                return_model = FALSE)

# This value is the BE/LE
be_lsd

par(mar = c(5, 5, 0.5, 0.5))
plot(diff(wl_detrend, 30)~diff(baro_detrend,30), transducer[seq(1, nrow(transducer), 10)],
     pch = 20, las = 1,
     cex = 0.05, col = '#000000',
     ylab = 'wl differences',
     xlab = 'baro differences',
     asp=1)
grid()

The least squares value of r round(be_ls, 2) is not very similar least squares differences method of r round(be_lsd, 2). Why might this be? What happens if we change the difference_space value?

# test 1 to 2880 lag_space
difference_spaces <- 2^(0:13)

be_lsds <- sapply(difference_spaces, function(x) {
  be_least_squares_diff(dat = transducer, 
                        dep = 'wl_detrend', 
                        ind = 'baro_detrend',
                        lag_space = x,
                        inverse = FALSE,
                        return_model = FALSE)
})

par(mar = c(5, 5, 0.5, 0.5))
plot(be_lsds~difference_spaces, log = 'x', type = 'o', pch = 20)
abline(h = be_ls, lty = 2)
text(x = 1, y = be_ls, label = 'Least squares solution', pos = 4)
grid()

So by using a different lag_space we can get agreement between the two methods. Let's see how other methods perform.

-@Clark1967

This function has the same arguments as be_least_squares_diff. This method will likely suffer from a similar problem as we saw with be_least_squares_diff. Let's use differences of 1 day as this appeared to give consistent results to the original least squares method.

be_c <- be_clark(dat = transducer, 
                 dep = 'wl_detrend', 
                 ind = 'baro_detrend',
                 lag_space = 720,
                 inverse = FALSE,
                 return_model = FALSE)

be_c

-@Rahi2013

The Rahi method is similar to the Clark method but has a magnitude check in addition to the sign check that Clark's method does. The arguments are the same for the current implementation except we do not have a return_model option.

be_r <- be_rahi(dat = transducer, 
                dep = 'wl_detrend', 
                ind = 'baro_detrend',
                lag_space = 720,
                inverse = FALSE)

be_r

Ratio method [@Gonthier2007]

This method calculates the ratio for each difference. Then filters them based on the magnitude of the barometric pressure change (quant), and finally summarizes the data using stat. We will calculate the one-day differences, select those having the largest barometric pressure changes (top 10 percentile), and take the median of the resulting ratios.

be_rat <- be_ratio(dat = transducer, 
                   dep = 'wl_detrend', 
                   ind = 'baro_detrend',
                   lag_space = 720,
                   quant = 0.9,
                   stat = median,
                   inverse = FALSE)
be_rat

High-low method

The high-low method compares water level changes to barometric pressure changes for different time periods. We specify the column name that holds the time and a time_group which is the size of the time interval in seconds. Here we calculate the LE every day for our dataset. This method will tend to overestimate the BE/LE.

be_hl <- be_high_low(dat = transducer, 
                     dep = 'wl_detrend', 
                     ind = 'baro_detrend',
                     time = 'datetime',
                     time_group = 86400*7)

head(be_hl)

-@Smith2013

I like this method for data where there is a complex trend. It uses informed judgement to fit the best BE/LE value. We supply different values of BE which are then used to remove the barometric signal from the data. The goal is to have a smooth curve (or a result that matches our expectations) after the barometric pressure is removed.

be_visual(transducer[seq(1, nrow(transducer), 120)],
          dep = 'wl_detrend',
          ind = 'baro_detrend',
          time = 'datetime',
          be_tests = seq(0.6, 0.9, 0.05),
          inverse = FALSE,
          subsample = TRUE) %>%
  partial_bundle()

What value do you estimate for LE? How certain are you? Why is there stucture in the residuals? One reason is that this data is from a semi-confined well with non-negligible well-bore storage and we have violated the assumption that there is no time lag in the response. A better technique for this well would have been to use barometric/loading response functions.

Acworth method

@Acworth2008

@Acworth2016

@Acworth2017

@Rau2018]

This method uses frequency domain values and includes the Earth tide response. The BE/LE is calculated by using ratio of the barometric pressure harmonic (S~2~) and the water level harmonic (S~2~) corrected for Earth tides. It is meant to be applied to confined wells. This value is also tied to the semi-diurnal frequency and care needs to be taken on it's interpretation in the same way the time based static barometric efficiency calculations are limited. See the response function vignette for a more complete view of the water level response.

be_a <- be_acworth(transducer, wl = 'wl', ba = 'baro', et = 'et', 
                   inverse = FALSE)

be_a

References

jkennel/waterlevel documentation built on Dec. 1, 2019, 6:24 p.m.