SEA Working Paper 00/12
Explaining groundwater hydrographs: separating atypical rainfall events from time trends
Ruhi FerdowsianA, David J. PannellB , Clare McCarronA, Arjen RyderA and Lisa CrossingA
AAgriculture
Western Australia, Albany WA 6330, Australia.
BFaculty of Agriculture, University of
Western Australia, Perth, WA 6907, Australia.
Corresponding author: Ruhi Ferdowsian: rferdowsian@agric.wa.gov.au
Abstract
By 1994, an estimated 1.8 million hectares of cleared land in Western Australia was affected by secondary dryland salinity to some extent. The area affected is likely to double in the coming 20 years. The cause of this salinity is excessive recharge under traditional agriculture, leading to rising groundwater levels. Monitoring changes in groundwater levels is helpful in indicating the degree of threat to agricultural land and public assets. Many researchers have studied groundwater level rises and attempted to explain them statistically.
We present an approach for statistically estimating trends in groundwater levels. The approach separates the effect of atypical rainfall events from the underlying time trend and the lag between rainfall and its impact on groundwater is explicitly represented. Rainfall is represented as an accumulation of deviations from average rainfall. Application of the approach is demonstrated using data from 49 bores in Jerramungup Shire, Western Australia. The approach provides high explanatory power, particularly for deeper bores.
Additional keywords: salinity, monitoring, sustainability indicators
Introduction
By 1994, an estimated 1.8 million hectares of cleared land in Western Australia was affected by secondary dryland salinity to some extent, representing 9.4 percent of agricultural land in the state (Ferdowsian et al. 1996). This area is likely to double in the coming 20 years (to around 3.3 million ha) and may double again before a new equilibrium is reached (Ferdowsian et al. 1996). The cause of this salinity is excessive recharge under traditional dryland (non-irrigated) agriculture, leading to rising levels of naturally saline groundwater (Wood 1924; Ghassemi et al. 1995). As water levels come close to the soil surface saline groundwater will discharge causing soil salinity and contaminating water resources. The rate of rise is not always consistent and increases during years with more than average rainfall.
Monitoring changes in groundwater levels in response to management practices is helpful in indicating the degree of threat faced, and the necessary timing and scale of preventative treatments. Many researchers have studied groundwater level rises and attempted to explain them statistically. A common approach is to fit a linear time trend to the data (e.g. Fig. 1) or to the annual minima when a strong seasonal response is evident. A potential problem with this approach is that rainfall during the period of observation may not be typical of long-term rainfall levels, so that observed rates of groundwater rise may not be relevant for future projections. Variations in rainfall between years within the sample period will also affect observed groundwater levels and fitted trends.
Fig. 1. Hydrograph with simple time trend (bore 7)
Shao et al. (1999) suggested an approach based on fitting different linear trends in different segments of the data over the duration of the records. Their approach has the advantage that it allows for differences in rainfall levels in different periods. It is also a useful method for classifying the hydrographs of a region into similar groups for further investigation.
However, there are limitations associated with their method:
Other approaches rely on analysis of time series data (e.g. Box et al. 1994; Larocque et al. 1998). These approaches are able to explain some of the seasonal variation in groundwater levels if regular and frequent monitoring has occurred. However, they do not explain the variation due to atypical rainfall events or atypical annual rainfall, which are apparent in most data series. They are also adversely affected by missing observations, which are particularly common, for example, in data sets collected by farmers.
In this paper we present and illustrate a new approach we have called HARTT (Hydrograph Analysis: Rainfall and Time Trends). The method can differentiate between the effect of rainfall fluctuations and the underlying trend of groundwater level over time. Rainfall is represented as an accumulation of deviations from average rainfall, and the lag between rainfall and its impact on groundwater is explicitly represented. In the following sections, we explain the approach in detail, and present an example of its application.
Data
The Jerramungup district is located on the south coast of Western Australia (118° 30¢ - 118° 30¢ E, 33° 15¢ -34° 30¢ S). The basement rocks in the north and central parts of the area are composed of medium and even-grained granites of Archaean origin, intruded by dolerite dykes. The weathering profiles in these areas are mostly shallow to moderately deep (<15 m), with deeply weathered (>15 m) profiles on the broad crests. The soils in the area are shallow (A horizon <0.25 m) duplex soils (Northcote 1979). The groundwater flow systems in the northern and central parts are of a local scale and the aquifers are separated by granite highs. Tertiary sediments with shallow to moderately deep (A horizon <0.5m) duplex soils cover the southern parts of the area. A complete Tertiary profile could be >50 m deep and consists of Pallinup Siltstone overlying the Werillup Formation. The Werillup Formation consists of dark-coloured clay and silt (lignite) overlying medium to coarse sand. The Pallinup Siltstone was deposited over the Werillup Formation in a marine environment. Closed depressions and swamps occur in areas where the siltstone lies over the Werillup Formation. This area is internally drained and has intermediate or regional scale groundwater flow systems.
The annual rainfall is 550 mm in the southwest but gradually reduces to 350 mm in the northeast. The pluviometric regime is characterised by long episodes of widespread and moderate rainfall but occasional short episodes of torrential and localised rainfalls may occur.
Agriculture Western Australia, a State government agency, has drilled more than 130 bores in this region (e.g. Henschke 1982; Martin 1992). Farmers have monitored the majority of these bores since 1990. Agriculture Western Australia is the custodian of the data and provides feedback to individual farmers and catchment groups.
Forty-nine bores were selected for analysis excluding any which had some treatment for reducing recharge (e.g. establishment of deep-rooted perennial plants). These bores had sufficient data (19 - 33 records) and were under traditional agriculture (continuous cropping or pasture/crop rotations) throughout the monitoring period. The sampling period varied between 7 and 10 years. The intended sampling interval was three months, but actual intervals varied between bores and were, in some cases, irregular, varying from three months to nine months. An important requirement for our analysis is that it be able to cope with irregularly spaced data and numerous missing values, as this is typical of data sets obtained from farmers’ monitoring.
Rainfall data were obtained from the Data Drill, a web site belonging to Queensland Department of Natural Resources (http://www.bom.gov.au/silo). Data Drill provides spatially interpolated data from daily rainfall records provided by the Australian Bureau of Meteorology from its network of recording stations.
Analytical method
Two forms of accumulative residual rainfall were used and compared:
The first was accumulative monthly residual rainfall (AMRR; mm):
(1)
where Mi,j is rainfall in month i (a sequential index of time since the start of the data set) which corresponds to the jth month of the year, Mj- is mean monthly rainfall for the jth month of the year, and t is months since the start of the data set.
Figure 2 illustrates the calculation of AMRR for an example bore. The AMRR variable is calculated as the difference between the other two variables shown. The AMRR variable tends to have relatively low within-year fluctuations because, in calculating AMRR, the fluctuations in actual rainfall tend to be offset by seasonal variation in average monthly rainfall.
Fig. 2. Calculation of AMRR (an illustrative example). AMRR is calculated from a data series beginning in 1957
The second was accumulative annual residual rainfall (AARR; mm):
(2)
where A- is mean annual rainfall. Because A- is a constant, the fluctuations in Mi are not moderated as they are for AMRR, so AARR has higher within-year fluctuations. For this reason, it was expected to be well correlated with data from bores with shallow groundwater levels (less than three metres deep) which typically have seasonally fluctuating water tables.
For both AMRR and AARR, construction of the variables was based on a data set commencing in 1957. This pre-dates the earliest recording of depth to groundwater by 33 years, allowing long lag effects of rainfall on groundwater to be detected, if they occur. Lags of up to several years were investigated.
The regression model was:
Deptht = k0 + k1 ´ AMRRt - L + k2 ´ t (3)
where Depth is depth of groundwater below the ground surface, t is months since observations commenced, L is length of time lag (in months) between rainfall and its impact on groundwater, and k0, k1 and k2 are parameters to be estimated. Parameter k0 is approximately equal to the initial depth to groundwater, k1 represents the impact of above- or below-average rainfall on the groundwater level, and k2 is the trend rate of groundwater rise or fall over time.
In 10 cases where it produced models with a higher R2 (most of which were shallow bores), AARR was substituted for AMRR in the regression model. The value of L was estimated separately for each bore by selecting the value that resulted in the highest R2 for the regression. Thus L does not necessarily represent the lag until either the first impact or the largest impact of rainfall on watertable depth, but the lag that produces the highest statistical correlation. In many cases, L is longer than the first detectable impact.
Results and discussion
Results for three individual bores are illustrated: deep (Fig. 3), moderate (Fig. 4) and shallow (Fig. 5). They show trends typical of other bores of similar depth. Results were similar, regardless of geology. For example, bores located in tertiary sand plains are marked in Table 1. Although geologically distinct from the other bores, they provided similar statistical results.
The deep bore (Fig. 3) which is typical of bores in broad crests near catchment divides, has a rapidly rising water table. The pattern of rise is extremely well captured by the estimated regression model (R2 = 0.97). The rise does not follow a simple linear trend, but the deviations from linearity are explained almost completely by the AMRR variable. The best fitting lag length for this bore is 17 months.
Fig. 3. Hydrograph for a deep bore (bore 7)
The moderately deep bore (Fig. 4), which is located in mid-slope, has a lower rate of rise (less than half that of Fig. 3). AMRR and time again explain groundwater level well, although in this case there is some unexplained variation (R2 = 0.80). The data deviates substantially from a smooth linear trend, but the main variations are captured well by the model. The best fitting lag is 2 months.
Fig. 4. Hydrograph for a medium-depth bore (bore 18)
Figure 5 is representative of shallow bores which are located in lower slopes and close to discharge sites. This bore has a slightly falling underlying trend. R2 is lower again (0.71) than the deeper bores, while again the broad deviations from linearity are captured by the rainfall variable (AARR in this case). The best fitting lag is zero. The correlation coefficients for most of the shallow bores (<3.0 m deep; Table 1) were low, indicating that other parameters are affecting groundwater level fluctuations. Vicinity to discharge zones, evaporation rate and hydraulic gradient could be amongst these parameters.
Fig. 5. Hydrograph for a shallow bore (bore 40)
As well as illustrating the quality of fit from the model, these three bores also demonstrate well the potential bias in estimation of groundwater trends that occurs in a simple regression model based only on time. Using the model of Eqn 3, the estimated underlying time trends of groundwater change are 0.0221, 0.0079 and –0.0033 m/month for the 3 bores, whereas the simpler model results in estimates of 0.0225, 0.0071 and –0.0046 m/month. The differences occur despite there being a long period of monitoring (>9 years). Of course, if data is available for a sufficiently long duration, the bias involved in using a simple linear time trend will be small. A strength of the HARTT method is that it improves trend estimates even where long data sets are not available.
Table 1 shows the regression results for all 49 bores in the data set. The p values shown for the two sets of parameter estimates are for t tests that the parameters are different to zero. As would be expected, the null hypotheses are less likely to be rejected in shallower bores. For both variables, this is partly because there is more unexplained variation in the data for shallow bores. For the time variable, this is reinforced by the tendency for parameter estimates to be lower for shallower bores. That is, for shallower bores, because the trend of groundwater change is closer to zero, it is less likely that we would reject a null hypothesis that it equals zero.
Table 1. Results of OLS estimations
It's a big table, so stored on a separate page. To view it, go to: http://www.general.uwa.edu.au/u/dpannell/dpap0012t.htm (184K)
Figures 6-10 reinforce and more fully describe a number of trends that have already been identified in the discussion so far. Fig. 6 shows the relationship between initial groundwater depth and the best fitting time lag between rainfall and water table depth. There is a clear trend for increasing lag with increasing depth, with between-bore variation also increasing with depth. The estimated lag lengths may perhaps be considered surprisingly high. We note that the lags shown are not the lags until the first impact of rainfall, but the lags that provide the best statistical correlation. This means that they are likely to be close to the maximum impact. The delay in this maximum impact is partly because of the slowness of groundwater movement downslope from land above the bore.
Fig. 6. Best-fitting lag length (L) as a function of initial groundwater depth
There is considerable variation in the estimated parameters for AMRR and AARR variables (Fig. 7). The model results indicated that groundwater level rose by between 0.5 and 6 mm per mm of above-average rainfall (Table 1).
Fig.7. Rainfall parameter estimate (K1) as a function of initial depth to groundwater
Rates of groundwater rise tend to be substantially greater for deep bores than for shallow bores (Fig. 8). This partly reflects the greater role of discharge in landscapes with shallow bores; this discharge is effective in removing water from the system and lowering the rate of groundwater rise. Indeed, for the group of bores with initial depths below 5 m, the average rate of groundwater rise is approximately zero.
Fig.8. Rate of groundwater rise (K2) as a function of initial depth to groundwater
Figure 9 reinforces the previous observation that the model fits groundwater data more precisely for deeper bores (although it should be noted that the number of deep bores is not high). Conveniently, it is more important to us to be able to predict groundwater rise for deeper bores. Indeed, farmers with shallow bores, where we have greatest difficulty in explaining data series and predicting trends, have little use for this information for these bores. Where the data is most needed, the quality of model estimates appears very high.
Bore 4 was unusual in this respect, in having a low R2 despite an initial depth to groundwater of 17.17 m (Table 1 and Figure 9). The same bore was an outlier in Figures 10 and 11. We observed that several water level measurements in this bore were exceptionally high. This could very likely be due to poor sealing of the bore allowing some of the perched water in the shallow soil profile to spill into the bore.
Fig.9. R2 as a function of initial depth to groundwater
The inclusion of the AMRR or AARR variable in the statistical model substantially increases its explanatory power (Fig. 10). The difference is lowest for very deep or very shallow bores, but can be dramatic for intermediate depths.
Fig.10. Comparisons of R2 for proposed method and simple linear time trend
Strengths and weaknesses of the approach
The method has a number of clear strengths: simplicity; separation of underlying time trend from unusual or atypical rainfall effects; improved estimation of time trends even where long data sets are not available; consistency of results with intuitive hydrological explanations; high degree of fit to data for bores with depths greater than 5 m; capacity to make high quality predictions of depth to groundwater for a period ahead equal to the length of the lag estimated in the model; and capacity to make less precise predictions even further ahead based on forecasts of rainfall for the period.
These strengths make the method a useful tool in groundwater resource management. However, some limitations can also be identified: (i) Without the addition of further variables to the model, the quality of explanation of data for shallow bores is not high. (ii) Trend results depend on calculated mean levels of rainfall. If means have changed during or after the period of rainfall data, the power of the model to forecast future groundwater levels is reduced. Changes in the within-year distribution of rainfall may also reduce predictive power. (iii) Depth changes over time so L and the coefficients of AMRR and t would actually change during the period of estimation. This would require a more sophisticated approach to estimation than we have undertaken. (iv) For time series data like this, serially correlated error terms are a risk which has not been allowed for in the estimation process.
The standard test for serial correlation (the Durbin-Watson test for first-order autoregressive errors – Durbin and Watson 1950) requires data to be available for consecutive periods. To allow this test, we redefined our time period as a quarter (three months). This introduces a slight error to some bores as the observations do not coincide precisely with quarters. Results of Durbin Watson tests for each bore are shown in Table 2, together with results for Ramsey’s (1969) RESET test for omitted variables and Cook and Weisberg’s (1983) test for heteroskedasticity.
Table 2. Diagnostic test results for OLS models
It's a big table, so stored on a separate page. To view it, go to: http://www.general.uwa.edu.au/u/dpannell/dpap0012t.htm (184K)
For several bores there are indications of omitted variables. Interestingly, these are more likely to be deeper bores, for which the quality of model fit is very high. Therefore the usual serious concern with omitted variables, bias in parameter estimates, is not so pressing.
Similarly, several bores fail the test for homoskedasticity. This is a less serious statistical problem, since parameter estimates remain unbiased, although inefficient (Judge et al. 1982). Procedures to estimate a weighted regression for these bores could be used if desired.
Unfortunately, most bores either fail the Durbin-Watson test outright, or fall in the uncertain range where failure is suspected. Statistically, when serial correlation is present, ordinary least squares regression (OLS) results in (a) parameter estimates that are unbiased but not efficient and (b) biased standard errors (Judge et al. 1982). That is, the estimates of rainfall impacts and time trends may be acceptable, but statistical inferences about them may be unreliable. For this reason, we explored use of standard transformations to attempt to overcome serial correlation.
Prais-Winsten (1954) proposed a transformed regression estimator to correct for first-order serially-correlated residuals (similar to the transformation of Cochrane and Orcutt 1949). The transformed regressions were conducted using Intercooled Stata version 6. Table 3 shows results from these transformed regressions for each bore. Overall, the success of the transformation in overcoming serial correlation was only moderate. We also note that in bores with missing observations, the transformations for serial correlation require the dropping of additional observations, further compounding the problem of missing information. Furthermore, even if there is bias in standard errors, it would have to be a bias of enormous magnitude to alter the conclusions from t tests in most of the bores, since few of the p values are at all close to the critical value of 0.05. Given these considerations, and the fact that some packages do not provide the transformation as an automatic option, it is worth considering further the possible error involved in using OLS.
Table 3. Results of PRAIS estimations
It's a big table, so stored on a separate page. To view it, go to: http://www.general.uwa.edu.au/u/dpannell/dpap0012t.htm (184K)
In the following comparisons we use quarterly data in each case, to allow us to focus on the effect of the estimation process itself, free of any distortions involved in converting data to quarterly. Figure 11 shows differences between the estimated AMRR or AARR parameter for OLS and the Prais-Winsten transformation. Given that the OLS parameter estimates are unbiased, there is no guarantee that the Prais-Winsten results are more accurate. Nevertheless, we are interested in whether the use of a method that is, in one sense, statistically preferred would substantially alter our parameter estimates. The figure shows that for bores deeper than 5 m, it probably would not. The exception is a single bore, which also corresponds to the outlier in Fig. 9.
Fig.11. Difference in ordinary least squares estimate of rainfall parameter (for AMRR or AARR) relative to PRAIS estimate
Figure 12 shows the corresponding graph for the time parameter. In this case, the use of OLS appears at least as robust. The Prais-Winsten transformation seems unlikely to alter time parameter estimates for bores deeper than 2 m.
Fig.12. Difference in ordinary least squares estimate of time trend relative to PRAIS estimate
Overall, where statistical inferences are to be drawn and differences between OLS estimated p values and critical values are not large, it may be judged necessary to account fully for serial correlation, even if it entails the loss of further information due to missing observations. However, where the primary goal is estimation of parameters rather than statistical inference, or where statistical inferences appear likely to be robust even in the presence of serial correlation, it would appear reasonable to rely on OLS.
Conclusion
We have presented an approach to statistical modeling of hydrographs that appears to have some considerable strength. The HARTT method is simple to apply with standard regression methods. It provides high quality fits to observed data in all but shallow bores (for which trend estimation is of less interest in any case). It allows the separation of atypical rainfall events from the underlying time trend. Results are highly consistent with hydrological expectations. In particular, we have examined a sample of 49 bores in Jerramungup, Western Australia and observed the following. Ground water depth is positively related to the best-fitting time lag between atypical rainfall and the groundwater level, the rate of groundwater rise, and the R2 for the fitted regression model. It is (weakly) negatively related to the parameter relating groundwater level to atypical rainfall. The method seems likely to have wide applicability in situations with rising groundwater trends.
Acknowledgements
We thank Michael Burton, David Dole and Greg Hertzler for statistical advice and assistance, Don McFarlane for support, the farmers of Jerramungup for their participation in the bore monitoring program and three anonymous referees for helpful comments and suggestions. David Pannell acknowledges financial support from Grains Research and Development Corporation.
References
Box GEP, Jenkins GM, Reinsel GC (1994) ‘Time series analysis: Forecasting and control.’ Vol.3 (Prentice-Hall: Englewood Cliffs, NJ)
Cochrane D, Orcutt GH (1949) Application of least squares regressions to relationships containing autocorrelated error terms. Journal of the American Statistical Association 44, 32-61.
Cook RD, Weisberg S (1983) Diagnostics for heteroscedasticity in regression. Biometrika 70, 1-10.
Durbin J, Watson GS (1950) Testing for serial correlation in least squares regression. Biometrika 37, 409-428.
Ferdowsian R, George R, Lewis F, McFarlane D, Short R, Speed R (1996) The extent of dryland salinity in Western Australia. In ‘Conference Proceedings, 4th National Conference and Workshop on the Productive Use and Rehabilitation of Saline Lands, Albany, Western Australia, 25-30 March 1996’. pp. 89-97. (Promaco Conventions: Perth, Western Australia)
Ghassemi F, Jakeman AJ, Nix HA (1995) ‘Salinisation of land and water resources: human causes, extent, management and case studies.’ (University of New South Wales Press: Sydney)
Henschke CJ (1982) Preliminary groundwater investigation in relation to soil salinity at Fitzgerald, W.A. Technical. Report No. 9, Division of Resource Management, Western Australian Department of Agriculture, Perth.
Judge GG, Hill RC, Griffiths W, Lütkepohl H, Lee TS (1982) ‘Introduction to the Theory and Practice of Econometrics.’ (Wiley: New York)
Larocque M, Mangin A, Razack M, Banton O (1998) Contribution of correlation and spectral analysis to the regional study of a karst aquifer. Journal of Hydrology 205, 217-231.
Martin S (1992) Groundwater investigations in the Jerramungup Shire. Technical Report No. 122, Division of Resource Management, Western Australian Department of Agriculture, Perth.
Northcote KH (1979) ‘A Factual Key for Recognition of Australian Soils.’ (Rellim: Adelaide).
Prais SJ, Winsten CB (1954) Trend Estimators and Serial Correlation. Cowles Commission Discussion Paper No 383, Chicago.
Ramsey JB (1969) Tests for specification errors in classical linear least squares regression analysis. Journal of the Royal Statistical Society, Series B 31, 350-371.
Shao Q, Campbell NA, Ferdowsian R, O'Connell D (1999) Analysing Trends in Groundwater Levels, A report prepared for Agriculture Western Australia as Part of a National Land and Water Resources Audit Project. Report Number: CMIS 99/37.
Wood WE (1924) Increase in salt in soil and streams following the destruction of native vegetation. Journal and Proceedings of the Royal Society of Western Australia 10, 35-47.
Citation: Ferdowsian, R., Pannell, D.J., McCarron, C., Ryder, A. and Crossing, L. (2000). Explaining groundwater hydrographs: separating atypical rainfall events from time trends, (SEA Working Paper 00/12). http://www.general.uwa.edu.au/u/dpannell/dpap0012.htm
Later revised version:
Ferdowsian, R., Pannell, D.J., McCaron, C., Ryder, A. and Crossing, L. (2001). Explaining groundwater hydrographs: Separating atypical rainfall events from time trends, Australian Journal of Soil Research 39: 861-875.
The SEA News index is at http://welcome.to/seanews