Groundwater is considered one of our most important natural resources. Since the early to mid-twentieth century, groundwater withdrawals have drastically increased due to changes in irrigation technology  . Population increases and industry demands have put additional pressures on groundwater resources. Recent studies estimate that groundwater extraction across the globe has tripled in the last 50 years and is increasing in many countries at an annual growth rate of between 1% - 2%    .
Climate variability brings another interesting set of challenges for groundwater management. Changes in either precipitation volume or distribution can impact recharge  . Regions which experience an increase in annual precipitation may see lower recharge due to a larger percentage of precipitation occurring as high intensity events  . Another challenge is associated with groundwater’s perceived higher resiliency to climate variability. Regions which are experiencing greater evapotranspiration will notice an immediate impact on surface waters. In either an attempt to protect surface waters or ensure a more reliable water source, the expectation is that many users will turn to groundwater as their source of water    .
Increasing our knowledge of aquifer volumes and how this translates into groundwater levels and trends is a significant component of improving policies and regulations protecting this resource. Quantifying the impact of multiple demands upon an aquifer in the wake of climate variability remains a difficult task. Additionally, there is a lack of long term continuous groundwater level and precipitation data (>50 years) combined with high spatial resolution. Most observational well data are within a region which has been affected by pumping, and often monitoring was started after irrigation.
Much of the research on groundwater levels and groundwater level trends over time seem to identify either declines in groundwater associated with increased pumping, or groundwater level declines associated with increased pumping but exacerbated by decreased recharge. Russo et al. (2014)  , in a summary paper of groundwater levels across the United States, broadly describes groundwater levels as having decreased in most of the United States. Combing both unconfined and confined aquifers, Russo et al. (2014)  found that in general increased pumping was associated with lower groundwater levels. Most regional studies also support broader patterns, especially in drier regions. In the Edwards Aquifer in Texas, Loaiciag (2003)  found that groundwater pumping had a much larger impact on groundwater level declines than changes in climate.
Ng et al. (2010)  used climate predictions to examine recharge in semi-arid regions of the US. Their study focuses on regions like the High Plains aquifer and found that while precipitation rates may fluctuate by −20% to +20%, groundwater recharge is much more variable, −75% to +35% due to its often episodic nature. While the positive recharge values were associated with wetter years, the extent of recharge declines suggests that groundwater recharge is an extremely sensitive variable and decreases in precipitation will result in an amplification of groundwater declines  .
The management of groundwater requires the ability to translate monitoring well data, which is represented as elevation, to its impacts on surface water and what can be sustainably pumped. The problem of predicting well levels or associated recharge is challenging for two main reasons, the first is that the factors which affect groundwater recharge and groundwater well levels are numerous and include the following factors: topography, vegetation, geology, climate and land use  . The second problem is related to the first; due to the numerous factors involved in groundwater prediction, there remains the substantial lack of publications and data available   . In the past ten years, as high resolution digital spatial data has become available, the studies on groundwater recharge increased but are limited to a region or too large an area for local management.
Our study is located in Central Wisconsin. The location is ideal to study the impact of irrigation on groundwater levels due to the high density of high capacity irrigation wells. The groundwater supply in Central Wisconsin is vital to domestic water demands as well as those of agriculture, industry and municipalities. For example, three counties in Central Wisconsin, Portage, Adams and Waushara, use 78 billion gallons of ground- water per year. Of the 78 billion gallons, approximately 87% or 67% billion gallons of groundwater is used for irrigation  . Soil type is the main reason for such a heavy dependence on the groundwater supply. The majority of soils in this region are highly permeable sands and gravels resulting from past glaciations. These sandy soils have a low water holding capacity, which stores little moisture for plants  . Sandy soils discouraged irrigated agriculture until improved technology developed in the 1950’s, created inexpensive irrigation systems  . Since the 1950’s, irrigation has become a dominant feature in Central Wisconsin, and may be a reason for groundwater related stresses such as declines in surface and groundwater levels.
The topography influences land use and the location of irrigated agriculture. Currently irrigated agriculture is concentrated on flat sandy areas. Irrigated land covers about 40% of the area of interest (Figure 1) and contains 70% of Wisconsin’s high capacity wells. Other land covers include non-irrigated agriculture, coniferous and deciduous forests, grassland, scrubland and wetlands.
In some regions of Central Wisconsin, groundwater related stresses are reflected in surface water declines. In 2005-2009, reaches of the Little Plover River, a groundwater fed stream in Plover, Wisconsin, intermittently dried up  . Long Lake, a groundwater fed lake located near Plainfield, Wisconsin (32 kilometers south of Plover), has also dried up  . The most highly stressed surface water resources occur in areas where there is a greater amount of irrigation.
Suggested reasons for declines in surface and groundwater levels are intensive groundwater pumping and drier weather. Precipitation records from the National Oceanic and Atmospheric Administration (NOAA), combined with the Palmer Drought Index and the Standard Precipitation Index, indicate that Central Wisconsin has received close to average annual rainfall during the study period (2005-2009)   . Despite near average precipitation totals, questions remain about the effects and interactions of precipitation on groundwater levels.
Figure 1. The topography and location of high capacity wells in the central sands region.
Many studies have been conducted throughout the United States that relate groundwater pumping to declines in surface waters or decreases of water levels in monitoring wells      . Only recent permitting practices in Wisconsin have required high capacity wells to log pumping rates and volumes, so historic pumping data are rarely available. Despite this lack of data, the consumption of groundwater and its effects on the Central Sand Plain’s surface waters and groundwater levels has been studied extensively. Weeks and Stangland (1971)  examined the development of present and future irrigation in the sand-plain area and its effects on stream flow and groundwater levels in the late 1960’s. Gotkowitz and Hart (2008)  looked at groundwater consumption and land use in Waukesha Wisconsin. Clancy et al. (2009)  examined groundwater use and its potential effects on the Little Plover River in Plover Wisconsin, and Kraft and Mechenich (2010)  studied groundwater pumping and its effects on groundwater, lake, and stream flow levels in the central sands of Wisconsin. The relationship between groundwater pumping and declines in surface and groundwater levels is well established, but the interaction between changes in climate, groundwater withdrawals, and the water table response are not as well understood  .
The direct measurement of the surface and groundwater response to pumping is presumably complicated by changes in precipitation which have occurred in some parts of Wisconsin. Increases in precipitation in the central part of the United States were noted by Lettenmaier et al. (1994)  and McCabe and Wolock (2002)  . More recently Juckem et al. (2008)  compared time periods 1941-1970 to 1971-2000 and found that wetter conditions have occurred in southwestern Wisconsin from 1971-2000. These wetter conditions were thought to be the result of a sudden shift in precipitation called a “step increase”. This step increase in precipitation may have masked the true impact of groundwater pumping pressures in some areas of Central Wisconsin  .
The goal of this study was to determine the impact of long term irrigation groundwater pumping in the presence of changing precipitation. Without pumping logs or observation well data that records the “pre-irrigation” record, we needed to identify the threshold year or years where a significant shift in groundwater levels occurred. This allowed us to establish irrigation trends which were complicated by changing precipitation. We used statistical methods to extract the changing precipitation from groundwater irrigation impacts. This study is unique not only because it provides a statistical method for establishing irrigation groundwater withdrawals in the presence of changing annual precipitation, but additionally, results from one of the test areas (Plover region) were compared to a calibrated and peer reviewed published groundwater model   .
The area of interest in central Wisconsin known as the central sands (Figure 1), is a loosely-defined region characterized by a thick (often >30 m) mantle of sandy materials overlying rocks of low permeability. Landforms are composed of glacial outwash plains and terminal moraine complexes associated with the Wisconsin Glaciation. The area is approximately 11,200 square kilometers and is 31% cultivated crops  The eastern boundary was delineated using ecoregions  and glacial deposits  and the western border is the Wisconsin River. This study focuses on the “headwater” or upland part of the central sands, east of wetlands or drained wetlands. The region contains more than 80 lakes (>5 ha) and over 1000 km of headwater streams and wetlands, which are well connected to shallow, unconfined, sand and gravel aquifers  .
2.1. Groundwater Data
Groundwater data from four U.S. Geological Survey (USGS) monitoring wells were used for this study (Table 1) (USGS, 2009)  . Well names were given based on the locale or quadrangle. The four monitoring wells were chosen based on two rationales:
Table 1. USGS monitoring wells used for data analyses. Plover 1 represents the original well number and plover 2 is the replacement number.
the length and consistency of available records (Table 2), and their location within the study area (Figure 2). The Wautoma and Amherst Junction wells are located in areas with a low density of high capacity wells and are thought to be minimally influenced by pumping. Hancock and Plover are located in areas with a high density of high capacity wells and are thought to be impacted by pumping.
USGS monitoring well data consisted of daily automated measurements for Amherst Junction, Hancock, and Wautoma and monthly field measurements for Plover. Depth to water was subtracted from benchmarked elevations to obtain the well water elevations. Average monthly elevations were used in statistical analyses.
Table 2. Available data for USGS monitoring wells used in this study.
Figure 2. The location of monitoring wells and municipalities with climate stations within the study region.
Data from the monitoring well near Plover were recorded under two different well numbers. Well number 442623089302701 was used prior to April 14, 2006 and was replaced by well number 442622089302901. These well measurements were combined and referenced to a common datum. Both well numbers are represented in Table 1.
2.2. Precipitation Data
Three types of annual precipitation data were used in this study (Figure 2; Table 3): 1) Five weather stations in central Wisconsin; 2) an interpolated annual data set for the Wautoma location; and 3) composite precipitation data for division 5  . Long term monthly precipitation data (≥50 years), from the cooperative observer (COOP) station network, were accessed online through the National Climate Data Center  . Yearly values were used in statistical analyses and were calculated from monthly observations. Missing monthly measurements were interpolated using a weighted average of the three closest COOP stations.
The interpolated annual precipitation data set for Wautoma was calculated using the Inverse Distance Weighting Method    . This method uses a weighted average of annual precipitation from the 12 closest COOP stations. Closer stations carry a heavier weight and have a greater influence on the average. All stations were within 50 miles and had sufficient records. Interpolations were not made when there were less than 12 stations contributing to the data.
In addition to annual precipitation records, the Standard Precipitation Index (SPI) was chosen as the main precipitation variable in multiple regression models because it best represents the lag between precipitation and the response in monitoring well water levels. Monthly growing season (May-September) SPI data for Wisconsin climate division 5 were obtained from NCDC. The SPI is a normalized index that quantifies precipitation deficits, can be calculated for any desired duration, and takes into account time scales in the analysis of wet and dry periods for water availability and use   .
The 24-month SPI (SPI24) was used for this study based on the work done by Mayer and Congdon (2008)  . They found that the 24-month SPI has the least influence in regression equations during normal precipitation periods (i.e. when SPI values are close
Table 3. COOP climate stations within the study region and periods of record for the division 5 data and the interpolated data set at Wautoma.
to zero). Other precipitation variables such as moving averages or lags will have less influence during dry conditions when precipitation values are close to zero, and more influence under wet conditions as precipitation values get larger  . In this study the SPI24 improved the explanation and prediction of groundwater fluctuations in multiple regression models and was able to represent the systematic response to wet and dry periods that occurs at monitoring wells in the study region.
2.3. Statistical Analyses
To determine a period within a monitoring well’s record where the statistical parameters (mean and standard deviation) do not change with time, we needed to establish a stationarity period.
Stationarity describes a process in which natural systems fluctuate within an unchanging range of variability  . When non-stationarity develops, it indicates that a shift has occurred between the relationships of hydrologic data within a region. Non- stationarity can be caused by changes in data collection methods or physical changes, such as a fluctuation in precipitation, or water diversion like groundwater pumping   .
Stationarity may be difficult to detect when unknown variables or multiple variables influence the system. To recognize these impacts a “covariate” also plays an important role in this study. A covariate is a statistical term used to identify an interaction which is not measured but is observed in the record    . A covariate may be binary and is often referred to as either a hidden, lurking or dummy variable.
Groundwater pumping and changes in precipitation were thought to be the two main covariates affecting groundwater levels in Central Wisconsin. Observations of pumping and changes in precipitation have no records associated with their impact on groundwater levels; therefore, a binary data set was developed for each covariate. For example, when pumping was thought not to have an effect on the groundwater record the data set was defined as “off”. When pumping potentially began to impact groundwater levels, the data set was defined as “on”. Because the covariates are disconnected from the continuous groundwater data, they may or may not actually represent groundwater pumping or changes in precipitation.
Kendall’s tau trend test established spatial differences in annual precipitation throughout the study area and established regions where the covariate for the step increase in precipitation occurred. The Mann-Whitney test, which calculated a difference in median values, was used to determine if a step increase occurred. Bivariate analysis indicated when changes showed up in the groundwater record. Multiple regression with the analysis of covariance (ANCOVA) was used to quantify, explain and predict the changes due to precipitation or pumping on groundwater levels. Findings were corroborated with other independent study’s findings        .
The main objective of the statistical analyses was to predict, quantify and explain changes in groundwater levels due to pumping. Without any pumping data, potential pumping could only be expressed as a covariate within multiple regression models (where the pumping effect is either “on” or “off”) and the date of the switch is determined by looking for a threshold year. In addition to the possible pumping influence on groundwater levels, there was the added influence from changes in precipitation. Studies in Southwestern Wisconsin showed a step increase in precipitation in the early 1970s which affected base flow   . The step increase in precipitation is a statistically significant shift in the mean value over a short period of time (one to two years). A shift in precipitation could potentially complicate the analyses by masking the impact thought to be due to pumping. Because there was no data set for the step increase in precipitation, it was also treated as a covariate or binary variable; however, from Juckem et al. (2008)  it was expected that the “switch” would occur sometime in the early 1970s.
2.4. Trend Analysis
Trends in precipitation were evaluated to determine if and where changes in precipitation occurred in the study area. In addition, an increasing trend in precipitation during the same time period that declines were measured in monitoring wells indicated potential pumping impacts. Kendall’s tau, a nonparametric statistical technique, has been regularly used to examine linear trends in precipitation    . This nonparametric test was used because precipitation data were found to have a non-normal distribution.
Trends were calculated for annual precipitation. The period of record used was 1955-2008 which was based on the shortest precipitation record at the Montello COOP climate station (Table 3). Annual trend analysis included Wautoma’s interpolated precipitation data, the five COOP climate stations, and the composite central division data. Calculations were made online using the Free Statistics and Forecasting Software website  .
2.5. Mann-Whitney Test
The Mann-Whitney test, a non-parametric version of the t-test, was used to determine the presence of the Juckem et al. (2008)  1970/1971 step increase in precipitation. This test was used for annual precipitation records because these data contained outliers that skewed the distribution. The period of record used to find the difference in median precipitation values was 1933-2008. This produced a similar number of data points on either side of the 1970, 1971 time break (n = 38 years) and included the most recent records. The composite division 5 records, four of the five COOP stations and the interpolated data for Wautoma had annual precipitation data for the 1933-2008 period. Montello’s record did not start until 1955 so a shorter period of record had to be used (1955-1986; n = 16). The Mann-Whitney test was calculated in Mini Tab (version 15).
2.6. Bivariate Analysis
A bivariate analysis tests for a difference in the means of two linearly correlated data sets  . The bivariate technique uses time series measurements and has commonly been used to evaluate climate data such as precipitation, evaporation and temperature    . In this study, bivariate analyses were used to evaluate changes in groundwater levels at monitoring wells. To meet the requirement of normal data, yearly average depth-to-water measurements from monitoring wells were used in the analysis.
The bivariate analysis was used to find the threshold year when the potential pumping covariate started to affect monitoring well levels. This was accomplished by examining non-stationarity, or a change in mean, between correlated test and control monitoring well records. The results from the bivariate technique determined the year, direction, and magnitude of the change in mean caused by non-stationarity  .
The bivariate analysis uses a regional stationary series which consists of multiple stations around a test station. The regional series is assumed to be independent and free of systematic change  . Due to the lack of available monitoring well records, multiple stations could not be used to develop a regional stationary series. Therefore, individual control locations (Amherst Junction or Wautoma) were considered the stationary regional series, which was similar to the methods of Kirono and Jones (2007)  .
Bivariate analysis was initially tested on the control monitoring wells, Amherst Junction and Wautoma, to develop stationary periods of record for each well. The time period from 1958-2008 was used based on the shorter data set at Amherst Junction. These stationary periods were developed so that the control locations could serve as the regional series in further analysis with test monitoring wells.
Once the stationary periods were established at the control sites, the bivariate test was used to determine the threshold year possibly caused by the pumping covariate at the test monitoring wells: Hancock and Plover. Control sites located within the closest proximity to the test sites were used as the stationary data set.
2.7. Multiple Regression with ANCOVA
Multiple regression using ANCOVA was the primary statistical technique used for this research. The main reason for using this technique was to see if independent variables (precipitation and the dummy variables for pumping and the step increase in precipitation) could explain and predict the dependent variable (monitoring well water elevation)  . Slope coefficients from the final model equations quantified the magnitude of the changes. ANCOVA describes the addition of a covariate or dummy/binary variable into the regression model.
Monthly monitoring well water elevations for the growing season (May through September) were used for multiple regression analyses. The growing season months were chosen to limit complexity due to snowpack infiltration rates and because most groundwater use occurs during the growing season.
To achieve parsimony in the model, we restricted the model to the three variables which provided the best results: the SPI24, the dummy variable for the potential step increase in precipitation and the dummy variable for the potential increased impact of groundwater pumping. The SPI24 was used in the models as a proxy for actual precipitation and the step increase covariate was used at locations where it was determined with previous tests that a precipitation increase had occurred. Regression tests were processed with PROC REG in SAS version 8.2.
3. Results and Discussion
3.1. Annual Precipitation Trends
Trends in annual precipitation were examined to determine if spatial differences existed in the study region. The time period 1955-2008 was used for all precipitation data sets (the five COOP climate stations, the interpolated data set at Wautoma, and the composite central division records) because data from one of the COOP climate stations (Montello) began in 1955.
Figure 2 shows the results of the spatial difference in the annual precipitation trends, where circles represent no trend and triangles represent increased trends (significant decreasing trends were not found). Three precipitation stations, Hancock, Montello, and Wautoma, in the southern part of the study area show increased trends in annual cumulative precipitation while stations in the northern part of the study area (Stevens Point, Waupaca, Wisconsin Rapids) and the composite division data show no trend (Table 4).
Precipitation varies from year to year and affects groundwater levels and hydrologic flow paths especially if annual totals have been above or below average for long periods   . Long term declines in precipitation would help to explain decreases in surface and groundwater levels, but increased precipitation though time may hide the impacts of pumping.
3.2. Step Increase in Precipitation
At locations where precipitation was significantly increasing, a step increase was suspected. The Mann-Whitney test was used to determine whether a step increase in precipitation occurred between 1970 and 1971 at six climate stations and for the region.
Table 4. P-values from the Kendall’s tau trends test and regression equations for annual precipitation at the six precipitation stations and for the composite central division data for 1955-2008. P-value < 0.05 indicate a significant trend and + indicates that the direction of the trend is positive.
* indicates significant p-value < 0.05. + indicates the direction of the trend.
Records from 1933-2008 were used when the data was available. This period of record contained an equal number of observations on either side of the 1970-1971 break (n = 38) and included current precipitation records.
Three climate stations (Hancock, Montello and Wautoma) and the regional precipitation data set showed a significant increase in median annual precipitation between 1970 and 1971 (Table 5). These three climate stations are located in the southern portion of the study area. The increase in median annual precipitation in the south during the past 38 years may have helped disguise the increased groundwater pumping. For this reason, regression equations for monitoring wells in the south at Hancock and Wautoma required a dummy variable for the change in precipitation while the monitoring well at Plover (in the north) did not. The dummy variable became active between 1972 and 1973 presumably to account for the groundwater lag which was found after running several regression models.
3.3. Bivariate Analysis
3.3.1. Control Monitoring Wells
The bivariate test detects the year, magnitude and direction of a systematic change in the mean between a test series and a second correlated stationary series. Control locations at Amherst Junction and Wautoma were considered the second correlated stationary series and the test series because both control locations were thought to be influenced by the covariate that represented the step increase in precipitation and the pumping covariate. For this reason, the bivariate test was initially used to identify a period of stationarity between the control monitoring wells.
The control monitoring well at Wautoma was thought to the least influenced by anthropogenic processes (i.e., pumping) due to the low density of irrigation wells. The bivariate test was calculated using Wautoma as the test series and Amherst Junction as the second stationary series for 1958-2008. Amherst Junction, a control well not greatly influenced by pumping, is located in the northern part of the study area where no step increase in precipitation occurred between 1970 and 1971.
Table 5. Results for changes in the median cumulative annual precipitation using the Mann- Whitney test for before and after 1970. P-value < 0.05 indicates a step increase in precipitation.
* indicates significant p-values < 0.05.
In Figure 3(a), a single discontinuity in the mean at Wautoma occurred in 1973, the year after the peak in the graph (1972). The dashed horizontal line in Figure 3 represents the 95% critical value. The peak, To, represents the maximum value of the difference (Ti) between the Wautoma and Amherst Junction data series. To occurs the year before the change in mean  , therefore non-stationarity was interpreted to occur after To (after 1973).
Non-stationarity that occurred in 1973 at the Wautoma monitoring well indicated that the increase in precipitation contributed to an increase in groundwater levels. The change in stationarity at Wautoma with respect to Amherst Junction occurs about the same time as a suspected step increase in precipitation for the area. Similar results were found by Lettenmaier et al. (1994)  when they used the bivariate test to determine that increases in stream base flow could be connected to increases in precipitation.
Figure 3. Graphs (a) (Top) and (b) (Bottom) of bivariate results for a change in mean at Wautoma monitoring well using Amherst junction as the stationary data set for two different time periods: 1958-2008 and 1972-2008. The dashed line represents the 95% critical value, Ti is the difference in the two data series being tested, and the vertical line is the peak year (To) which occurred one year before the change in mean. This discontinuity in mean is associated with the step increase in precipitation between 1970 and 1971. Statistics (Ti) below this line indicate no change in mean, establishing a stationary period between 1972-2008.
A stationary period from 1972-2008 was established at Wautoma, which included the peak (1972), but excluded the years prior (1958-1971). Figure 3(b) illustrates that with the years prior to 1972 excluded, Ti does not reach the critical value, which suggests a change in mean did not occur at Wautoma for the new time period (1972-2008). The results for the new stationary period implied that the only change to the monitoring well at Wautoma had to do with the step increase in precipitation which occurred around 1970.
To find a stationary period at the Amherst Junction control well, Amherst Junction was used as the test series and Wautoma was the second correlated stationary series. In the Wautoma record it was found that non-stationarity was associated with the step increase in precipitation. Amherst Junction is located in the northern part of the study area where there was a small or no step increase in precipitation (Waupaca p-value = 0.066 in Table 5). Therefore, the entire Wautoma record (1958-2008) was used to establish a stationary period at Amherst Junction.
Two shifts in mean occurred in the Amherst Junction record. Although the bivariate test was designed to detect a single change in mean, it can be sensitive to multiple changes with the largest shift identified as the primary break and the smaller shift identified as the secondary break  . The first change in mean occurred in 2000 after the peak in 1999. The bivariate test was reevaluated without 2000-2008 and identified another change in mean which occurred during in the early 1960’s. Three peaks were found in the years 1961, 1962, and 1965. Multiple peaks meant that from 1962-1966, there were multiple adjustments and responses occurring at the monitoring well. To identify the longest possible stationary period at Amherst Junction, data prior to the first change in mean (1962) were left out. When 1958-1961 were excluded from the bivariate calculation, a stationary period from 1962-1999 was established.
Using the bivariate test, two stationary periods were found in the records of the control monitoring wells. At Wautoma the stationary period was from 1972-2008 and at Amherst Junction it was from 1962-1999. The stationary periods at the control wells were important to establish because they were used to detect a change in mean at the test monitoring wells. Non-stationarity at the test monitoring wells indicated the year when a threshold was reached where groundwater pumping was suspected to have a measureable impact on groundwater levels.
3.3.2. Test Monitoring Wells
To find the threshold year for the potential pumping covariate at test monitoring wells, bivariate analysis was calculated using the control monitoring wells as the stationary data set. The bivariate analysis determined the year, magnitude and direction of non-stationary periods which may have been caused by pumping. Control wells closest in distance to test wells were compared due to similar precipitation patterns.
The Wautoma control well was used to find the implied pumping covariate threshold at Hancock, because both were influenced by the step increase in precipitation and because of proximity. Wautoma’s stationary period (1972-2008) did not include the time period before 1970. After 1971, the step increase in precipitation had already occurred. Therefore, the effect of precipitation did not influence the determination of the pumping threshold year at Hancock.
The outcome of the bivariate test at Hancock indicated that although 1994 was the peak year, the results plateau from 1994-1998. Potter (1981)  mentions that while a plateau does not create a clear estimate of the exact year of change, it demonstrates that the bivariate test is sensitive to any change taking place. The plateau may coincide with a ramp up of groundwater pumping between 1994 and 1998. The last peak in the plateau (1998) was selected as To which resulted in a non-stationarity, or a second stationary period after 1999. Multiple regression models corroborated 1999 as the possible pumping threshold year. The bivariate analysis results for Hancock indicate that an increase in the magnitude of groundwater pumping during the mid to late 1990’s may have caused declines in groundwater levels.
The Plover test well was initially compared to the Amherst Junction control well (1962-1999) because Amherst Junction was the closest control location. The bivariate results for the comparison show two peaks in the Plover record between 1962 and 1999. The first peak, in 1973, signified non-stationarity at Plover starting in 1974. A second, higher peak occurred in 1986 indicating an additional non-stationary period starting in 1987. In a previous study, using double mass curves, Clancy et al. (2009)  determined that groundwater declines became noticeable in the Little Plover River around 1973. The first peak was similar to results from Clancy et al. (2009)  . Therefore, 1973 was chosen as the first potential pumping threshold year instead of 1986, even though 1986 represented a slightly higher peak.
The multiple peaks indicated a possible second non-stationary period, so the Plover test well was compared with the Wautoma control well. The second comparison to the Wautoma control well was used to determine whether additional declines in groundwater levels occurred later in the record. The Wautoma stationary period (1972-2008) started close to the first identified threshold year at Plover (1973). Therefore, the earlier non-stationary peak found when Plover was compared to Amherst Junction did not influence the comparison of Plover to Wautoma.
The bivariate test using Wautoma as the correlated stationary series identified an additional non-stationarity period at Plover between 1972 and 2008. The year 1999 was identified as the second possible pumping threshold because it represented the end of the plateau period or a period of continuous change.
3.4. Multiple Regression with ANCOVA
Multiple regression equations were developed for the three monitoring wells. Regression models used the growing season (May-September) water elevations from 1960-2008. Variables used in the equations included the 24-month Standard Precipitation Index (SPI24), and the dummy variables for the step increase in precipitation (STP) and groundwater pumping (P1 and P2). Equations for each monitoring well are:
Hancock = 0.35∙SPI24 + 0.48∙STP + −0.97∙P1 + 325.74 R2 = 0.77 Equation (1)
Plover = 0.33∙SPI24 + −0.39∙P1 + −0.89∙P2 + 330.49 R2 = 0.69 Equation (2)
Wautoma = 0.20∙SPI24 + 0.36∙STP + 264.78 R2 = 0.63 Equation (3)
Amherst Junction = 0.42∙SPI24 + 1.40∙STPC-0.92∙PC1 + 338.73 R2 = 0.57 Equation (4)
Monitoring well locations are water elevations (m), SPI24 is the standard precipitation index for 24 months, STP is the dummy variable for the step increase in precipitation (m), P1 is the first dummy variable for pumping (m), and P2 is a second dummy variable for pumping (m). All regression models were significant with p-values < 0.05.
To better illustrate monitoring well responses to the SPI24, STP, P1 and P2, the Multiple Regression and ANCOVA section is broken down into subsections based on each monitoring well.
3.4.1. Case Study 1: Hancock Monitoring Well
At the Hancock monitoring well the increase in groundwater levels potentially due to the STP after 1973 was 0.48 m. When the difference between observed and predicted groundwater levels was observed in 1999, P1 was added to the model and the calculated decline was 0.97 m (Equation (1)).
Figure 4(a) shows the predicted and observed groundwater levels at Hancock for the growing season. The predicted water elevations include the STP dummy variable but no pumping variable. The vertical lines represent dates the STP variable was activated and
Figure 4. Graphs (a) (Top) and (b) (Bottom) of observed and predicted multiple regression results at the Hancock monitoring well for the growing season (May-September) 1960-2008. Graph (a) includes the STPC after 1972 and graph (b) includes PC1 which began to affect monitoring well levels in 1999.
represents irrigation pumping impacts in the record. Figure 4(b) includes P1 and shows the predicted results of Equation (1). The P1 variable modified the suggested precipitation response to give a more accurate prediction of the observed monitoring well water levels.
We consider the withdrawal of groundwater to be the main reason for declines at the Hancock monitoring well during the end of the record because annual and seasonal precipitation totals were higher than during the previous time periods. Although pumping was developed before 1999, the STP may have reduced the measureable impact pumping had on groundwater levels. Due to lack of pumping records and because increases in precipitation may be masking pumping at Hancock, it is difficult to calculate the full effect of pumping on this location.
3.4.2. Case Study 2: Plover Monitoring Well
We assumed that the Plover monitoring well was not influenced by the increase in precipitation that occurred at Hancock in the early 1970s because it is located in the northern part of the study region. As there was no real increase in precipitation and a steady increase in irrigated agriculture, pumping was thought to influence groundwater levels at an earlier date.
Before 1974, precipitation was thought to be the dominant influence of water levels in the Plover monitoring well, although pumping developed in the area prior to the 1970’s. As shown in Figure 5(a), the first response in the Plover regression model, which may have been associated with pumping (P1), occurred after 1974 and resulted in a water level decline of 0.39 meters. Around 1998 the model was no longer able to predict water levels. The additional groundwater decline (P2) after 1998 was 0.89 m (Equation (2)) and is shown in Figure 5(b). The net decline at the Plover monitoring well for both periods was 1.28 m. At Hancock the decline was only slightly lower (0.97 m) and both wells seemed to respond to a possible increase in pumping at around the same time (1999).
3.4.3. Case Study 3: Wautoma Monitoring Well
The regression model for the monitoring well at Wautoma included only the SPI24 and the STP. The well is located in the southern part of the study region where the increase in precipitation occurred and it is in an area with very little irrigated agriculture suggesting little pumping influence.
Simplifying the Wautoma regression model by using only SPI24, the model adjusted its predicted values to the time period after 1972. Figure 6(a) shows that predicted water elevations before 1972 were lower. When the STP was included the predicted water elevations adjusted downward before 1972 (Figure 6(b)). The STP added 0.36 m to the modeled monitoring well water levels (Equation (3)). This was similar to the modeled response from the step increase in precipitation at Hancock (STP = 0.48 m).
The Wautoma monitoring well water elevations show little fluctuation and the regression model predicts more sharp peaks and troughs than what actually existed in the data. This lack of response to seasonal fluctuations throughout the record may be due
Figure 5. Graphs (a) (Top) and (b) (Bottom) of observed and predicted multiple regression results at the Plover monitoring well for the growing season (May-September) 1960-2008. Graph (a) includes PC1 which occurred after 1973 and graph (b) includes PC2 added after 1998.
to Wautoma’s position in the groundwater flow system which is lower than the two other monitoring wells. Additionally, the Wautoma well does not respond quickly to precipitation events. This is illustrated with the low modeled SPI24 slope coefficient of 0.20 (Equation (3)). The rate of change is lowest among the three monitoring well equations. The slow response might explain the model’s attempt to include nonexistent peaks and troughs.
3.4.4. Case Study 4: Amherst Junction Monitoring Well
The Amherst Junction monitoring well responded to two shifts in the record between 1960 and 2008. At the end of the record, after 1999, water elevations in the monitoring well declined by 0.92 m (Equation (4)). This was similar to declines found during the same time period at Plover (0.89 m) and Hancock (0.97 m). Although Amherst Junction is thought to be influenced by groundwater pumping, the monitoring well is located in an area with fewer high capacity wells and therefore the influence from pumping was thought not to be as great. Pumping may have contributed to the decrease in water levels, but the magnitude of the response might have been caused by less precipitation. The reference to this precipitation shift is lower summer precipitation totals which occurred at the Waupaca climate station from 2000-2008.
Figure 6. Graphs A (Top) and B (Bottom) of observed and predicted multiple regression results at the Wautoma monitoring well for the growing season (May-September) 1960-2008. Graph (a) shows the response to the SPI24 before the STP was added. Graph (b) includes the STP added after 1972.
At the beginning of the record there was a positive shift in the Amherst Junction water elevations that did not correlate with the step increase in precipitation which occurred after 1972. The increase in water elevations after 1962 was 1.40 m. In Equation (4), this variable is referred to as STPC, but it does not correspond to the same time period used for the STPC at either Hancock or Wautoma.
A possible explanation for the shifts in the Amherst Junction monitoring well water levels before 1962 and after 1999 could partially be due to the well’s location. The Amherst Junction well is located on the shores of Lake Emily in western Portage County. During three different occasions in the well’s record values were measured above the land surface. The location of the monitoring well could also explain the large fluctuations in water elevations through the record and during the growing season.
Figure 7. Graphs (a) (Top), (b) (Middle) and (c) (Bottom) of observed and predicted multiple regression results at the Amherst Junction monitoring well for the growing season (May-Sep- tember) 1960-2008. Graph (a) shows just the SPI24, graph (b) includes the water level decline that occurred after 1999 and graph (c) contains the increased water levels after 1962.
the SPI24 precipitation variable, clearly shows the breaks in the record before 1962 and after 1998. Figure 7(b) and Figure 7(c) include the addition of the two other variables (STPC and PC1). The predicted values from the model are improved, although it seems as though the variations in Amherst Junction’s record are too large for the model to accurately explain.
3.5. Multiple Regression Summary
The difference between predicted regression values and observed water elevations were thought to result from pumping or changes in precipitation patterns. An earlier modeled pumping signature at Hancock may have been masked by the step increase in precipitation. Predicted regression values at Plover may show an earlier pumping signature because there was no increase in precipitation. The cumulative effect of the step increase at Hancock was thought to be similar to that found at Wautoma, even though Wautoma showed little fluctuation in groundwater levels. Table 6 quantifies changes in groundwater levels at each monitoring well.
Surface and groundwater levels have declined in some regions of the study area possibly due to groundwater withdrawals. Groundwater levels have also increased in other regions due to a suggested step increase in precipitation. Annual and seasonal trends revealed that precipitation increases occurred in the southern part of the study area and were generally associated with summer rainfall. In the northern part of the study region, no significant trend in precipitation was detected so a spatial difference between the northern and southern part of the study area had to be taken into account when examining the potential pumping/precipitation interaction. The Mann-Whitney test confirmed trend tests and identified those locations where a step increase in precipitation occurred between 1970 and 1971.
In order to extract the influence of precipitation changes to suspected pumping on groundwater well levels, binary covariate variables were required in multiple regression equations. These variables represented the potential impacts of pumping and precipitation and revealed declines in groundwater levels. An increase in precipitation added an average of 0.42 m to groundwater levels at Hancock and Wautoma. At Hancock, the increase in groundwater levels was thought to mask the effects of pumping earlier in the record, making the quantification of groundwater declines before 1999 difficult. The net decline in groundwater levels at the Plover monitoring well was 1.28 m.
The conclusions confirm the hypothesis for this study. Increases in precipitation have changed monitoring well levels by increasing groundwater levels in some regions of the study area. Groundwater levels have declined in other regions of the study area despite increases in precipitation. The use of multiple statistical approaches and the
Table 6. Results from multiple regression models which quantify increases and declines in monitoring well water elevations (m) possibly due to pumping or the step increase in precipitation. Dates represent when the dummy variables were added to the regression model.
STP: Dummy variable for the step increase in precipitation. P1: Dummy variable for the first pumping signature. P2: Dummy variable for the second pumping signature.
corroboration with recent studies by Clancy et al. (2009)  and Kraft and Mechenich (2010)  give a strong inference that there is limit to the sustainability of surface and groundwater systems.
This research was funded by the Wisconsin Department of Natural Resources (project NMI00000247).