Modeling of GRACE-Derived Groundwater Information in the Colorado River Basin

Groundwater depletion has been one of the major challenges in recent years. Analysis of groundwater levels can be beneficial for groundwater management. The National Aeronautics and Space Administration’s twin satellite, Gravity Recovery and Climate Experiment (GRACE), serves in monitoring terrestrial water storage. Increasing freshwater demand amidst recent drought (2000–2014) posed a significant groundwater level decline within the Colorado River Basin (CRB). In the current study, a non-parametric technique was utilized to analyze historical groundwater variability. Additionally, a stochastic Autoregressive Integrated Moving Average (ARIMA) model was developed and tested to forecast the GRACE-derived groundwater anomalies within the CRB. The ARIMA model was trained with the GRACE data from January 2003 to December of 2013 and validated with GRACE data from January 2014 to December of 2016. Groundwater anomaly from January 2017 to December of 2019 was forecasted with the tested model. Autocorrelation and partial autocorrelation plots were drawn to identify and construct the seasonal ARIMA models. ARIMA order for each grid was evaluated based on Akaike’s and Bayesian information criterion. The error analysis showed the reasonable numerical accuracy of selected seasonal ARIMA models. The proposed models can be used to forecast groundwater variability for sustainable groundwater planning and management.


Introduction
The portion of water stored in pore spaces of soil and rock beneath the Earth's surface is known as groundwater.Groundwater is a major and dependable large source of freshwater for domestic, agricultural, and industrial users in all climatic regions of the world [1,2].It also plays a central role to maintain the balance in the ecosystem [3,4].The areas with highly variable rainfall fully depend on groundwater due to lack of adequate alternative resources of fresh water.The over-exploitation of groundwater is a threat to the sustainability of ecosystems, water supply, and economic and social developments [5].More energy is needed to pump out the dwindling groundwater levels, which triggers higher energy consumption [6].Ground subsidence may also occur for the excessive groundwater depletions [7,8].Long-term sustainable groundwater from river basins has become a major concern among scientists.It is difficult to quantify the rate of natural renewal of groundwater to subsidize the water table depletion required for supporting the ecosystem [9].A proper understanding of groundwater storage (GWS) variation at spatial and temporal scale helps to understand the hydrologic cycle and its effect on global climate change.Therefore, consistent monitoring of the groundwater levels and storage may help to develop better plans for maintaining balanced ecosystems, sustainable economic development and encountering the water stress [10].
Routine monitoring of groundwater from regional to continental scale with the networks of monitoring well is tedious, time-consuming, and expensive.The inconsistency in data collection at the temporal scale, scarce of evenly distributed monitoring wells, the complicated subsurface physical properties, and recharging processes make the groundwater head estimation more complicated.Sometimes to access the hydrological data is restricted by the law enforcement agency.Since 1970s scientists are working on satellite-based remote-sensing systems for measuring hydrological component; groundwater is the latest addition on that [11].Avoiding aforementioned challenges, National Aeronautics and Space Administration (NASA)'s twin Gravity Recovery and Climate Experiment (GRACE) mission provides the first opportunity to directly measure groundwater changes from space through the accurate estimation of Earth's gravity field variations over time and space [12][13][14].Unlike traditional remote sensors utilizing electromagnetic emissions, GRACE uses K-Band microwave to estimate inter-satellite drift due to the gravitational force, the effect of the redistribution of earth mass component including cryosphere, hydrosphere, ocean, atmosphere, and land surface within the accuracy of a micrometer.The GRACE data is processed with numerical models to truncate the effect of atmospheric and oceanic contributions.Therefore, GRACE-observed mass changes mostly reflect terrestrial water storage (TWS), vertical integration of groundwater, soil moisture (SM), surface water (SW), snow water, and vegetation water measurements.Removing the change in surface water storage from GRACE gravity measurements yields an estimated change in groundwater [6,[15][16][17][18].GRACE provides groundwater level variability relative to the user-defined time baseline instead of an optimal value of groundwater level from the surface or any datum.GRACE-derived GWS demonstrated higher potential to represent the in-situ groundwater-level [18][19][20][21].Many studies have already quantified groundwater variability from GRACE gravity measurements [22][23][24].
Analyzing the historical variability of the GRACE-derived GWS will help to understand future groundwater conditions.Trend analysis is one of the popular and most accepted methods of analyzing the variability of any data and is also popular while analyzing the changes in a hydrologic variable [25][26][27].Hydro-meteorological data are often non-normally distributed where mean does not represent the central tendencies as effectively as the median.In such cases, non-parametric tests are often useful over parametric tests to analyze the time series data.Mann-Kendall (MK) a commonly used non-parametric statistical test to detect the trend in hydro-meteorological data [25,26,28].MK test was introduced by Kendall, [29] and Mann, [30] and is thus named after them.It is the rank-based trend detection analysis which gives the statistical significance of the trend if they are present in the time series.Further, the magnitude of the trend is often evaluated with the Thiel-Sen approach in most of the time series analysis [26,28,31].
Modeling and predicting of groundwater level fluctuation provides valuable information regarding groundwater declination, trend and allowable limit of exploitation.Forecasting the potential change in groundwater levels can help to mitigate the issues associated with water management [32].Accurate prediction of groundwater head helps planners to augment urban, rural, and industrial freshwater demand.Groundwater level forecasting also helps to comprehend the dynamics between factors that affect groundwater tables.An accurate and reliable groundwater forecasting helps to develop integrated management of groundwater and SW [5].The comprehensive river basin management policies require a continuous water level assessment, modeling, and forecasting.Over the past years, the conceptual and physically-based numerical models are widely used in groundwater simulation and quantification.These models are able to replicate the physical properties of groundwater dynamics.However, these modeling techniques have practical limitations due to the absence of long-term high-quality data and the complex structures of aquifers [33][34][35][36].In such cases, statistics-based time series analysis can be used as a suitable alternative [33,37].The time series analysis is able to deal with the data that are collected at a regular interval, like groundwater measurements.Time series analysis is widely used in groundwater resources management [38][39][40][41].
The prediction process in time series analysis highly depends on previous observations.Among various statistical approaches, a stochastic process is commonly used to capture the trend of forecasted data considering uncertainty.Autoregressive Integrated Moving Average (ARIMA) model is one of the most popular models that follow the behavior of the stochastic process [42].ARIMA and autoregressive moving average model (ARMA) is often utilized to understand and forecast different time series data.Unlike ARMA, ARIMA model is beneficial while working with the non-stationary time series data.ARIMA model considers the underlying correlation and patterns in the lagged data to forecast future data.ARIMA is the combination of differencing, auto-regression (AR), and moving average (MA) technique.The ARIMA deals with the fewer coefficients, which is the foremost benefit of using this model.ARIMA is the mathematical approach for predicting the future scenario considering the changes in trends and serial correlation among the previous observations.The use of the ARIMA model to forecast the hydrological time series data are well documented.In hydro-climatological research, the ARIMA model was used for predicting the mean monthly streamflow [43,44], rainfall [39,45], long-term runoff [46], river water quality and discharge [47,48], and drought [49].
The need for sustainable decisions through groundwater assessment and prediction has motivated the current study.The inconsistency of in-situ groundwater data and the sparseness of such data at higher spatiotemporal scale motivated utilizing the GRACE-derived groundwater anomaly in the current study.Additionally, deterministic forecasting models can mimic the spatial variability of groundwater while they undergo parametric uncertainties resulting in poor forecasts.Thus, to subside such shortcomings, GRACE-derived groundwater anomaly was coupled with ARIMA model for forecasting groundwater anomaly.The study was tested in one of the drought-affected area-the Colorado River Basin (CRB).The Colorado River is governed by the 'Law of River' and suffered severe water stress during recent climate change and increasing water demand.Addition to surface water, decline the reductions in groundwater storage was observed.The study tests the following research questions: (1) Is GRACE-derived data applicable to analyze the groundwater variability in the region undergoing drought like CRB? (2) What are the historic spatiotemporal variations in the groundwater in the CRB? (3) Can a stochastic ARIMA model coupled with GRACE data forecast future groundwater variability at the selected spatiotemporal scale of the CRB?
The aforementioned research questions are addressed by utilizing GRACE-derived groundwater changes as such data provides a longer length of data at uniform spatiotemporal scales.Before using the GRACE data for further study, it was tested against the in-situ observations to verify its reliability.The spatiotemporal variability of the historic GRACE-derived groundwater anomaly was tested using non-parametric MK test, and the magnitude of the trend was analyzed with Thiel-Sen's approach.Then, the stochastic ARIMA model was developed by determining its parameters to accommodate the seasonality and autocorrelation in the groundwater.Finally, the ARIMA model was tested to simulate the historical and future groundwater variability in the CRB.The documented research reveals that significant efforts had been made in the past to forecast groundwater conditions especially in the regions with higher water stress.According to documented literature and the author's knowledge, this is the first time ARIMA is integrated with GRACE data to predict future groundwater anomalies.ARIMA model showed promising result while simulating historical groundwater records, which can be utilized for groundwater management.Associating the proposed model with other deterministic models may enhance the confidence of both approaches.This novel approach may help water managers in making sustainable groundwater policies.

Study Area
The study area comprising of CRB including major rivers is illustrated in Figure 1.CRB has an area of 657,000 km 2 which sources water supply to approximately 40 million people of seven different states (Colorado, Utah, Nevada, Arizona, California, Wyoming, and New Mexico), and Mexico [50].CRB is divided into the lower CRB and upper CRB, where 8.6 million and 1 million people reside respectively [51].Streams mostly originate in the upper part of the basin approximately ninefold as compared to the lower part.This basin provides water to a total of 22,000 km 2 of irrigated land [50].Population growth and increased irrigation in the CRB has boosted water consumption.Increasing water demands and severe drought since 2000 have posed the declination of water levels in the CRB's largest reservoirs-Lake Powell and Lake Mead.To offset paucity of surface water supplies in drought, water users have switched to using groundwater to a greater extent.But the irregularities of groundwater reserves in CRB raises questions: how much water is already consumed, and how long can it be sustained?Castle et al. (2014) [52] showed that the CRB experienced a loss of 41 million acre-feet groundwater from 2003 to 2014, 77 percent of total freshwater decrement.Groundwater pumping for irrigation triggered the losses in water storage [53].Groundwater is being depleted at a much faster rate than previously thought indicating that the CRB is going to experience a groundwater decline in the coming years.

Data Sources
The data sets used in the analysis comprises of GRACE and Global Land Data Assimilation System (GLDAS) data [54].Section 3.1 provides details about the GRACE data followed by Section 3.2 which discusses how different hydrologic components of water balance equation were obtained from GLDAS data.Section 3.2 also presents the assimilation of GRACE and GLDAS data to calculate groundwater storage anomaly.

Data Sources
The data sets used in the analysis comprises of GRACE and Global Land Data Assimilation System (GLDAS) data [54].Section 3.1 provides details about the GRACE data followed by Section 3.2 which discusses how different hydrologic components of water balance equation were obtained from GLDAS data.Section 3.2 also presents the assimilation of GRACE and GLDAS data to calculate groundwater storage anomaly.

GRACE Data
GRACE is twin satellite launched by the NASA, and the German Aerospace Centre combined in March 2002 for tracking down the mass redistribution of the earth by monitoring changes in gravitational force [14].As GRACE record cumulative signals both gravitational and non-gravitational effects, atmospheric and oceanic effect need to be isolated to quantify the change in hydrological mass [55,56].GRACE provides monthly data at both 1 • × 1 • [57] and 0.5 • × 0.5 • [58] spatial resolution.To reduce noise and measurement errors, several filtering processes such as a destriping filter, a 200 km wide Gaussian averaging filter, and the low-pass spectral filter are applied [13].The filtering process for getting smoother data weakens true geophysical signals.A set of scaling factors is considered with the data set to restore the signal [59].The present study uses GRACE Level-3 RL 05 anomaly data, with temporal and spatial resolution of one month and of 1 • × 1 • respectively, obtained from the Center for Space Research at the University of Austin/Texas (CSR), NASA Jet Propulsion Laboratory (JPL) and the German Research Centre for Geosciences (GFZ) to calculate GWS variations.While computing the equivalent water height, the average of CSR, JPL, and GFZ was taken to reduce the noise [21].The time series of variability in gravity field for each cell from January 2003 to July 2016 were used in this study.This dataset has gaps of several months.The imputeTS [60], an R package, was used to take care of missing values as the ARIMA model cannot deal with the time series with missing data [61].The TWS comprises GWS, SM, SW, snow water equivalent (SWE), and canopy water storage (CWS).The scaling factor was used to restore the information of the GRACE-derived product, ignoring the effect of the SW component [62].The change in TWS was calculated using the water balance equation as shown in Equation (1).

Global Land Data Assimilation System (GLDAS)
GLDAS is a joint project of NASA, the National Oceanic and Atmospheric Administration, and the National Centers for Environmental Prediction which simulate hydrologic components by integrating ground-based and satellite observations at a higher temporal and spatial resolution [54].GLDAS provides hydrological data using four Land Surface Models (LSMs), i.e., the Community Land Model (CLM) [63], Variable Infiltration Capacity (VIC) [64], Noah [65], and Mosaic [66].To measure the changes in GWS from GRACE TWS, the values of SM, SWE, and CWS need to be deducted from the land surface model of GLDAS.The average of four LSMs models (CLM, VIC, Noah, and Mosaic) was used to calculate the anomaly of SM, SWE, and CWS to minimize any biases or errors [20].All LSMs with a similar spatial and temporal resolution of GRACE was used in this study.The GLDAS derived products (SM, SWE, and CWS) were also converted to the same anomaly of GRACE.
SM/SWE/CWS anomaly at time t, was evaluated as shown in Equation (2).
where P represents SM/SWE/CWS and P 2004−2009 represents corresponding average P from January 2004 to December 2009.Finally, groundwater storage anomalies (GWSA) was calculated utilizing Equation (3) shown below.
The in-situ groundwater head was obtained from monitoring stations, which was a measure of relative height from North American Vertical Datum of 1988.Groundwater level change from different monitoring wells was factored by respective specific yield to get GWS.Then, GWSA from the observed value was calculated using the same equation used for GLDAS data.This GWSA had been evaluated in terms of change in groundwater level represented hereafter as ∆GWL.Adjusted ∆GWL from in-situ was compared with GRACE-derived ∆GWS.

Methodology
First, Section 4.1 discusses the evaluation of historical variations in groundwater storage.Followed by Section 4.2 which presents the method of forecasting groundwater anomaly with the ARIMA model.

Evaluation of Trend in the Groundwater Data
MK-a rank-based test-was utilized to detect the presence of the trend in the GRACE-derived ∆GWS data.The test also revealed the statistical significance of the trend if any in each gridded data.More details on the MK test can be obtained from Mann [30] and Kendall [29].The current study only considered the trend if they were detected at 5% statistical significance based on the MK test.The trend magnitude was calculated utilizing the Theil-Sen approach [31].

ARIMA Model
The ARIMA model framework is shown in Figure 2. First, the groundwater data was procured as presented in Section 3. Figure 2 illustrates the steps involved in developing an ARIMA model corresponding to the GRACE-derived groundwater anomaly of one grid which was repeated for the groundwater analysis of entire CRB.As shown in Figure 2, first the groundwater anomaly was derived from the GRACE data which was tested for the presence of stationarity.ARIMA model parameters were derived from the stationary data and if the data was not stationary it was derived from the differenced data.Once, accurate model parameters were obtained, the historical inputs were then utilized to forecast the groundwater anomaly.The equations involved in building an ARIMA model is discussed later in this section.
Before evaluating the model parameters, stationarity in the data was evaluated.Differencing was done to make the data stationary, and the process was accomplished for each grid as shown in Figure 1.The non-seasonal ARIMA model with first-order differencing is expressed as Equation ( 4) shown below.
In general, an ARIMA model is represented as ARIMA (p, d, q) where, p, d, and q are the order of autoregression, integration (differencing), and moving average order, respectively.The more simplified ARIMA model is mathematically expressed as shown in Equation ( 5) where B is the backward shift operator.
Most of the hydrological data follow a seasonal pattern, like monthly, quarterly etc.The autocorrelation function (ACF) and partial autocorrelation function (PACF) test were performed to determine the seasonality in the dataset.A time series data, attribute with the seasonality, is converted to stationary in nature with the help of seasonal differencing.A seasonal ARIMA model is classified as an ARIMA (p,d,q) (P,D,Q) m model, where P = order of the seasonal autoregression, D = order of seasonal differencing, Q = order of moving average and m = the number of periods in one season.
considered the trend if they were detected at 5% statistical significance based on the MK test.The trend magnitude was calculated utilizing the Theil-Sen approach [31].

ARIMA Model:
The ARIMA model framework is shown in Figure 2. First, the groundwater data was procured as presented in Section 3. Figure 2 illustrates the steps involved in developing an ARIMA model corresponding to the GRACE-derived groundwater anomaly of one grid which was repeated for the groundwater analysis of entire CRB.As shown in Figure 2, first the groundwater anomaly was derived from the GRACE data which was tested for the presence of stationarity.ARIMA model parameters were derived from the stationary data and if the data was not stationary it was derived from the differenced data.Once, accurate model parameters were obtained, the historical inputs were then utilized to forecast the groundwater anomaly.The equations involved in building an ARIMA model is discussed later in this section.Before evaluating the model parameters, stationarity in the data was evaluated.Differencing was done to make the data stationary, and the process was accomplished for each grid as shown in For the backward shift operator, B, a seasonal ARIMA (p,d,q) (P,D,Q) m can be expressed as follows: where, ϕ p (B) is Auto regressive operator, θ q (B) is moving average operator, Φ p (B m ) is Seasonal autoregressive operator, and Θ q (B m ) = Seasonal moving average operator.ARIMA models require stationary time series data for modeling.The stationarity of the input time series data needs to be identified first.If the data series has any non-stationarity component, the non-stationarity should be removed before utilizing the data as ARIMA model inputs.The ACF and PACF plots are used to recognize if a time series exhibits non-stationary or not.ACF shows the correlations between time series and its own lag.PACF also depicts the correlations between the time series and its own lag along with neglecting the effect members stay in between them.The two parallel lines, with the x-axis, in the plots remarks 95% confidence intervals (CI) for the calculated ACF and PACF.Besides an additive seasonal decomposition is performed to visualize the seasonal patterns of time series data.
An additive seasonal decomposition can be written as where at any period t, seasonal component is represented by S t , T t is the trend component, and R t is the remainder component Differencing is usually used to stationarize the nonstationary data time series.The lowest differencing order for which time series maintain a well-defined mean value and ACF plot decays rapidly to zero, either from above or below, is the appropriate order for differencing.ACF and PACF generally give an idea of tentative order for moving average and autoregressive of ARIMA model.Different order of combinations close to the potential model was generated.
Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) were utilized to choose the appropriate ARIMA model.Equations ( 8) and ( 9) below present the expressions to evaluate AIC and BIC.
where the sum of squared errors represented by SSE is calculated as in the statistical model, N is the observations number and ε represents the white noise.The minimum values of these AIC and BIC criteria signify better model performance [67].
The errors in the model are often accessed with the residuals.Residuals are random in nature which cannot be explained by mathematical interpretations such as distribution with zero mean and constant variance.They are uncorrelated and normally distributed.ACF and PACF of the residuals are standard evaluation criteria for ARIMA model evaluation.For normally distributed residuals, all the spikes in ACF and partial ACF should range within the threshold values [68].The fitted model with minimum ACF and PACF for each gridded data were used to forecast the groundwater anomalies.3) was used to estimate the monthly distributed groundwater storage anomalies based on GRACE using the above-mentioned hydrological components.GRACE-derived ∆GWS was correlated with the average of in-situ measurements in both upper and lower CRB with the R 2 value of 0.62 and 0.67 respectively.The R 2 values testify good conformity between satellite and ground-based measurement as for such analysis the R 2 between 0.55 and 0.75 is ranked as good correlation [21].Like previous studies, the larger study area leads to getting better co-relation among GRACE and in-situ observations.

Historical Variations in Groundwater
From Figure 3a, it can be observed that the groundwater storage in the upper basin declined sharply during 2013.This can be attributed to the lowest recorded snowfall in Rocky Mountain and extreme drought scenario in 2012 [52].In lower CRB as observed in Figure 3b, groundwater storage followed a decreasing trend from 2004 to 2012 before it experienced a significant depletion from 2012 to mid-2014.On the other hand, the rise in groundwater storage was visible for a shorter period between mid-2009 to mid-2010).Moderately wetted weather lessened the demand for surface water supplies and groundwater recharged during this period.The consistency in groundwater head, in both upper and lower basin, after January 2015 shows strong evidence of recovery from the drought starting from 2000 up to 2014. Figure 3c shows the spatial variation of the trend at the statistical significance of p ≤ 0.05.The magnitude of the trend was calculated as the Theil-Sen's slope is represented in color on the map.The trend magnitude was computed for each grid independently as shown in the figure.This magnitude represents the change in the GRACE-derived ΔGWS in centimeters per year (cm/year).The negative trend magnitude suggests the decline in ΔGWS over time while the positive trend magnitude signifies the increase in ΔGWS over time.As evident from Figure 3c, most of the grid (80 grids out of 94 grids) had shown the significant decline in ΔGWS at the statistical significance of p ≤ 0.05.The maximum negative trend magnitude observed in the CRB was 2.34 cm/year.Thirteen grids did not show any trend at the statistical significance of p ≤ 0.05, and only one grid of lower CRB showed the increase in ΔGWS at the statistical significance of p ≤ 0.05 with the trend magnitude of 0.34 cm/year.As observed in Figure 3c, although entire CRB has experienced a decline in ΔGWS, upper CRB has undergone more decline in the ΔGWS as compared to the ΔGWS of lower CRB.This might be attributed to the recent management and conservation policies in the lower CRB and varying recharge of upper and lower CRB [69].The decline in ΔGWS deems the assessment of the future groundwater variability in the region.

ARIMA Model Results:
The ARIMA model project future ΔGWS data based on the past trend in GRACE-derived ΔGWS.Selecting the order of ARIMA model is an iterative process which is discussed in the later part of the paper.The absence of constant mean, variance, and autocorrelation over time make the time series data non-stationary.Generally, the non-stationarity appears in two aspects, i.e., changes in the parameters and the feature of the driving cause.Both are difficult to be identified without the prior assumed way the non-stationarity happens.Thus, the stationarity of time series data should be tested Figure 3c shows the spatial variation of the trend at the statistical significance of p ≤ 0.05.The magnitude of the trend was calculated as the Theil-Sen's slope is represented in color on the map.The trend magnitude was computed for each grid independently as shown in the figure.This magnitude represents the change in the GRACE-derived ∆GWS in centimeters per year (cm/year).The negative trend magnitude suggests the decline in ∆GWS over time while the positive trend magnitude signifies the increase in ∆GWS over time.As evident from Figure 3c, most of the grid (80 grids out of 94 grids) had shown the significant decline in ∆GWS at the statistical significance of p ≤ 0.05.The maximum negative trend magnitude observed in the CRB was 2.34 cm/year.Thirteen grids did not show any trend at the statistical significance of p ≤ 0.05, and only one grid of lower CRB showed the increase in ∆GWS at the statistical significance of p ≤ 0.05 with the trend magnitude of 0.34 cm/year.As observed in Figure 3c, although entire CRB has experienced a decline in ∆GWS, upper CRB has undergone more decline in the ∆GWS as compared to the ∆GWS of lower CRB.This might be attributed to the recent management and conservation policies in the lower CRB and varying recharge of upper and lower CRB [69].The decline in ∆GWS deems the assessment of the future groundwater variability in the region.

ARIMA Model Results
The ARIMA model project future ∆GWS data based on the past trend in GRACE-derived ∆GWS.Selecting the order of ARIMA model is an iterative process which is discussed in the later part of the paper.The absence of constant mean, variance, and autocorrelation over time make the time series data non-stationary.Generally, the non-stationarity appears in two aspects, i.e., changes in the parameters and the feature of the driving cause.Both are difficult to be identified without the prior assumed way the non-stationarity happens.Thus, the stationarity of time series data should be tested before performing any further analysis.First, the intermediate sample results of ARIMA forecast corresponding to the grid with centroid 110.5 • W and 31.5 • N referred hereafter as sample grid is presented.After presenting the results corresponding to the sample grid the results pertaining to all grids of CRB are presented.The ACF and PACF correlograms are drawn in order to determine the stationarity in monthly ∆GWS data as shown in Figure 4a  If the ACF and PACF promptly converge to zero with the increase in the time lag, the time series is considered to be stationary.As seen in Figure 4a, the spikes in the ACF plot are over the threshold values, and their net values are decreasing with time.This prominently indicates the presence of nonstationarity in the groundwater.Besides, the spikes in the PACF plot don't shut off immediately.In the PACF plot of Figure 4b, several numbers of spikes are above the threshold.This indicates the presence of seasonality in the data.After the occurrence of a trend or non-homogeneity, the time series changes its regime and no longer remains stationary.It is hard to model such non-stationary time series; thus the stationary part of the time series should be extracted before creating any model.ARIMA is associated with the difficulty to model and forecast the time series which are subjected to any changes in the regime.Therefore, 12-month differencing was applied to make it stationary by removing the seasonal influence.The seasonally differenced time-series data are shown in Figure 4c.
Figure 5a,b shows the ACF and PACF plots for seasonally differenced time series corresponding to the sample grid.These plots provide a tentative order of the ARIMA model.The ARIMA model with the tentative order is termed hereafter as the sample model.From Figure 5a, tentative model order cannot be obtained as the higher number of spikes are above the threshold in the ACF plot.There is no significant spike in the PACF as shown in Figure 5b plot for non-seasonal lags, suggesting a possible MA model of zero order.Whereas, for the seasonal component, there is a significant spike at lag 12 signifying the first order seasonal AR corresponding to the sample grid.From this, a preliminary seasonal ARIMA (p, 0, q) (1, 1, 0)12-the potential model was selected to simulate the inter-annual variability GRACE-derived ΔGWS from January of 2003 to December of 2014.In total twelve models corresponding to the sample grid, with several combinations of orders were assumed, If the ACF and PACF promptly converge to zero with the increase in the time lag, the time series is considered to be stationary.As seen in Figure 4a, the spikes in the ACF plot are over the threshold values, and their net values are decreasing with time.This prominently indicates the presence of non-stationarity in the groundwater.Besides, the spikes in the PACF plot don't shut off immediately.In the PACF plot of Figure 4b, several numbers of spikes are above the threshold.This indicates the presence of seasonality in the data.After the occurrence of a trend or non-homogeneity, the time series changes its regime and no longer remains stationary.It is hard to model such non-stationary time series; thus the stationary part of the time series should be extracted before creating any model.ARIMA is associated with the difficulty to model and forecast the time series which are subjected to any changes in the regime.Therefore, 12-month differencing was applied to make it stationary by removing the seasonal influence.The seasonally differenced time-series data are shown in Figure 4c.
Figure 5a,b shows the ACF and PACF plots for seasonally differenced time series corresponding to the sample grid.These plots provide a tentative order of the ARIMA model.The ARIMA model with the tentative order is termed hereafter as the sample model.From Figure 5a, tentative model order cannot be obtained as the higher number of spikes are above the threshold in the ACF plot.There is no significant spike in the PACF as shown in Figure 5b plot for non-seasonal lags, suggesting a possible MA model of zero order.Whereas, for the seasonal component, there is a significant spike at lag 12 signifying the first order seasonal AR corresponding to the sample grid.From this, a preliminary seasonal ARIMA (p, 0, q) (1, 1, 0) 12 -the potential model was selected to simulate the inter-annual variability GRACE-derived ∆GWS from January of 2003 to December of 2014.In total twelve models corresponding to the sample grid, with several combinations of orders were assumed, close to the tentative order to estimate ARIMA model parameters.The comparison of AIC and BIC of different models are summarized in Table 1.As shown in the table, all the model has the same orders of differencing, as the model was selected minimizing AIC and BIC [70].Among all models, ARIMA (1,0,1) (1,1,1) 12 showed the minimum AIC and BIC.Residual analysis of a fitted model depicts the robustness of the model.For a well-fitted model, the residuals should behave as white noise (random).The ACF and PACF of ARIMA (1,0,1) (1,1,1) 12 model residuals for the sample grid are plotted in Figure 5c,d respectively.As shown in Figure 5c,d, spikes of both ACF and PACF, lies within 95% confidence level, indicating the randomness in residuals.Besides, as all the lags have coefficients below the threshold, the absence of autocorrelation between them is confirmed.1.As shown in the table, all the model has the same orders of differencing, as the model was selected minimizing AIC and BIC [70].Among all models, ARIMA (1,0,1) (1,1,1)12 showed the minimum AIC and BIC.Residual analysis of a fitted model depicts the robustness of the model.For a well-fitted model, the residuals should behave as white noise (random).The ACF and PACF of ARIMA (1,0,1) (1,1,1)12 model residuals for the sample grid are plotted in Figure 5c,d respectively.As shown in Figure 5c,d, spikes of both ACF and PACF, lies within 95% confidence level, indicating the randomness in residuals.Besides, as all the lags have coefficients below the threshold, the absence of autocorrelation between them is confirmed.6, the ARIMA model forecasts were close to the GRACE-derived ∆GWS during the testing period (January 2015-December 2016) signifying good model skills.The model also captured the seasonality and the trend of ∆GWS series remarkably well during the testing period as shown in Figure 6.In the training period, the simulated ∆GWS followed the decreasing trend of 0.43 cm/year which was close to the decreasing trend of GRACE-derived ∆GWS evaluated as 0.47 cm/year.On the other hand, during the testing period, both ARIMA results and GRACE-derived ∆GWS showed the recharge in groundwater.Tillman et al. [71] used Coupled Model Intercomparison Project Phase 5 climate projections and supported the idea of recharge in groundwater storage of the CRB.It was observed that the model captured most of the peaks with some minor exceptions.The incompetence of the model to capture some of the peaks may be attributed to the presence of randomness component in the stochastic model.The selected seasonal ARIMA model generates the groundwater storage variability with higher accuracy for the remaining time series.The higher value of Nash-Sutcliffe efficiency (NSE) (0.89) and comparatively lower Root Mean Square Error (RMSE) (1.04) during the training period and 0.72, 1.5 during testing period respectively, indicating higher conformity between satellite-derived and simulated groundwater anomaly for the sample grid.This also suggests robustness of the seasonal ARIMA model that offer anticipated consistency and accuracy.Finally, the tested ARIMA model was utilized to forecast future groundwater anomaly of the sample grid for 36 months.
ARIMA (2,0,1) (  The model also captured the seasonality and the trend of ΔGWS series remarkably well during the testing period as shown in Figure 6.In the training period, the simulated ΔGWS followed the decreasing trend of 0.43 cm/year which was close to the decreasing trend of GRACE-derived ΔGWS evaluated as 0.47 cm/year.On the other hand, during the testing period, both ARIMA results and GRACE-derived ΔGWS showed the recharge in groundwater.Tillman et al. [71] used Coupled Model Intercomparison Project Phase 5 climate projections and supported the idea of recharge in groundwater storage of the CRB.It was observed that the model captured most of the peaks with some minor exceptions.The incompetence of the model to capture some of the peaks may be attributed to the presence of randomness component in the stochastic model.The selected seasonal ARIMA model generates the groundwater storage variability with higher accuracy for the remaining time series.The higher value of Nash-Sutcliffe efficiency (NSE) (0.89) and comparatively lower Root Mean Square Error (RMSE) (1.04) during the training period and 0.72, 1.5 during testing period respectively, indicating higher conformity between satellite-derived and simulated groundwater anomaly for the sample grid.This also suggests robustness of the seasonal ARIMA model that offer anticipated consistency and accuracy.Finally, the tested ARIMA model was utilized to forecast future groundwater anomaly of the sample grid for 36 months.Similar to the sample grid, the analysis was performed for each grid corresponding to the GRACE data within the CRB. Figure 7 shows the distribution of ARIMA model orders for each grid within the CRB. Figure 7i,ii represents the non-seasonal and seasonal components of the ARIMA model respectively.In Figure 7, columns (a), (b), and (c) represents autoregressive order, differencing order, and moving average order respectively.It can be seen that the moving average and autoregressive order both in seasonal and no-seasonal part varies between 0 and 2. The differencing order for the seasonal part was found 1 in almost all the grids.Non-seasonal differencing order was adopted as 1 if the seasonal differencing order was 0 and 0 if the seasonal differencing order was 1.
within the CRB. Figure 7i,ii represents the non-seasonal and seasonal components of the ARIMA model respectively.In Figure 7, columns (a), (b), and (c) represents autoregressive order, differencing order, and moving average order respectively.It can be seen that the moving average and autoregressive order both in seasonal and no-seasonal part varies between 0 and 2. The differencing order for the seasonal part was found 1 in almost all the grids.Non-seasonal differencing order was adopted as 1 if the seasonal differencing order was 0 and 0 if the seasonal differencing order was 1.The ARIMA model was then established for each grid with the order as summarized in Figure 7.The resulting model results were promising to forecast future groundwater variability.Figure 8a shows the monthly distribution of ARIMA forecasts during the historical period of January 2003 to December 2016 while Figure 8b shows the historical records based on GRACE-derived data for the same period and within the CRB.Under detailed observation, it can be seen that there are 164 box plots in each figure representing the monthly inter-quartile range of ΔGWS variability spatially within the CRB.The whiskers on each box plots represent the corresponding 5th and 95th percentiles.While comparing Figure 8a,b, it can be observed that the corresponding box plots are similar and the monthly variations along the abscissa are also consistent.This validates the robustness of the ARIMA for simulating the GRACE-derived ΔGWS.The boxplots in Figure 8c,d   The ARIMA model was then established for each grid with the order as summarized in Figure 7.The resulting model results were promising to forecast future groundwater variability.Figure 8a shows the monthly distribution of ARIMA forecasts during the historical period of January 2003 to December 2016 while Figure 8b shows the historical records based on GRACE-derived data for the same period and within the CRB.Under detailed observation, it can be seen that there are 164 box plots in each figure representing the monthly inter-quartile range of ∆GWS variability spatially within the CRB.The whiskers on each box plots represent the corresponding 5th and 95th percentiles.While comparing Figure 8a,b, it can be observed that the corresponding box plots are similar and the monthly variations along the abscissa are also consistent.This validates the robustness of the ARIMA for simulating the GRACE-derived ∆GWS.The boxplots in Figure 8c,d  The median NSE during the training period was observed to be greater than 0.9 and that during the testing period was approximately 0.7.This shows the model holds very good standings as suggested by [72].As the ARIMA model showed skillful results during the historical period from the robustness of the ARIMA model is established by NSE during both training and testing period which varies between 0.5 and 1.The median NSE during the training period was observed to be greater than 0.9 and that during the testing period was approximately 0.7.This shows the model holds very good standings as suggested by [72].As the ARIMA model showed skillful results during the historical period from January 2003 to December 2016, the model was utilized to forecast future groundwater anomaly from January 2017 to December 2019.Figure 9a,b summarizes the forecasted monthly GRACE-derived ΔGWS in centimeters (cm) for upper CRB and lower CRB respectively.Each box plot shows the distribution of a monthly forecasted ΔGWS of entire grids independently for both upper and lower CRB.The boxes represent the corresponding 25th, and 75th percentile ΔGWS and the red lines represent median ΔGWS.Similarly, the whiskers represent the 5th and 95th percentiles of the corresponding ΔGWS.If the boxes lie in the positive ordinate, it shows the month will have higher groundwater storage as compared to the corresponding historical months.From Figure 9a,b, future groundwater in lower CRB is expected to decline more as compared to upper CRB as the distribution of ΔGWS is more negative for lower CRB represented by the boxes with minimal ordinates for lower CRB.It is the reflection of the ARIMA's property of making forecasts based on the historical variations as, during the recent historical drought period, groundwater within arid and drought-affected regions CRB was withdrawn to supplement the scanty surface water.Further, the ARIMA model was able to capture the seasonality in the groundwater anomalies.From Figure 9a, it was observed ΔGWS showed positive anomaly for most of the grids in upper CRB in the months like March, April, and May as the median ΔGWS was Figure 9a,b summarizes the forecasted monthly GRACE-derived ∆GWS in centimeters (cm) for upper CRB and lower CRB respectively.Each box plot shows the distribution of a monthly forecasted ∆GWS of entire grids independently for both upper and lower CRB.The boxes represent the corresponding 25th, and 75th percentile ∆GWS and the red lines represent median ∆GWS.Similarly, the whiskers represent the 5th and 95th percentiles of the corresponding ∆GWS.If the boxes lie in the positive ordinate, it shows the month will have higher groundwater storage as compared to the corresponding historical months.From Figure 9a,b, future groundwater in lower CRB is expected to decline more as compared to upper CRB as the distribution of ∆GWS is more negative for lower CRB represented by the boxes with minimal ordinates for lower CRB.It is the reflection of the ARIMA's property of making forecasts based on the historical variations as, during the recent historical drought period, groundwater within arid and drought-affected regions CRB was withdrawn to supplement the scanty surface water.Further, the ARIMA model was able to capture the seasonality in the groundwater anomalies.From Figure 9a, it was observed ∆GWS showed positive anomaly for most of the grids in upper CRB in the months like March, April, and May as the median ∆GWS was generally positive.This can be attributed to the groundwater recharge from early snowmelts as a result of changing the climate, in the snow-covered regions of upper CRB.Most of the grids especially in lower CRB showed negative anomalies for all months.This shows the groundwater condition of the CRB is unlikely to improve in the near future as compared to the historical records.This can be attributed to higher evaporation and lower precipitation during dry seasons resulting from climate change.Further, the decline in surface water availability also causes stress on groundwater.Once the groundwater storage is declined resulting from excessive withdrawal, it causes soil consolidation reducing the recharge rate of groundwater.These factors combined hinders the improvement in groundwater conditions.As the forecasted groundwater anomalies are mostly negative during August and September, it is more likely that groundwater storage will decline during these months of near future.
generally positive.This can be attributed to the groundwater recharge from early snowmelts as a result of changing the climate, in the snow-covered regions of upper CRB.Most of the grids especially in lower CRB showed negative anomalies for all months.This shows the groundwater condition of the CRB is unlikely to improve in the near future as compared to the historical records.This can be attributed to higher evaporation and lower precipitation during dry seasons resulting from climate change.Further, the decline in surface water availability also causes stress on groundwater.Once the groundwater storage is declined resulting from excessive withdrawal, it causes soil consolidation reducing the recharge rate of groundwater.These factors combined hinders the improvement in groundwater conditions.As the forecasted groundwater anomalies are mostly negative during August and September, it is more likely that groundwater storage will decline during these months of near future.

Conclusions
The study presents the stochastic behavior of GRACE-derived groundwater time series data in CRB from January 2004 to December 2016.An ARIMA model is developed using the training data from January 2004 to December 2014 to forecast the monthly time series data from January 2015 to December 2016.A total of eleven seasonal ARIMA models with a different order of combination are tested to find the best fit model.The key findings of the current study are summarized below: (1) GRACE-derived groundwater anomaly being well correlated with in-situ groundwater data established the applicability of GRACE data in analyzing past and future groundwater analysis in the CRB.
(2) The GRACE-derived change in monthly groundwater storage showed strong seasonality.Seasonal differencing was promising while making forecasts with non-stationary groundwater data.
(3) The ACF and PACF plots of differencing data series were beneficial to estimate the tentative order of the ARIMA model.
(4) ARIMA models order obtained based on the evaluation criteria (AIC and BIC) were found skillful.The residual analysis reinforced the idea of selection of the fitted model order.
(5) ARIMA estimates of groundwater storage anomalies fit reasonably well with the observed values as supported by RMSE and NSE skills during the historical training and testing periods.(6) The ARIMA forecasts indicated the increase in March, April and May groundwater storage within the major number of grids of upper CRB.This can be attributed to the early snowmelts in the region during these months as a result of climate change.

Conclusions
The study presents the stochastic behavior of GRACE-derived groundwater time series data in CRB from January 2004 to December 2016.An ARIMA model is developed using the training data from January 2004 to December 2014 to forecast the monthly time series data from January 2015 to December 2016.A total of eleven seasonal ARIMA models with a different order of combination are tested to find the best fit model.The key findings of the current study are summarized below: (1) GRACE-derived groundwater anomaly being well correlated with in-situ groundwater data established the applicability of GRACE data in analyzing past and future groundwater analysis in the CRB.(2) The GRACE-derived change in monthly groundwater storage showed strong seasonality.Seasonal differencing was promising while making forecasts with non-stationary groundwater data.(3) The ACF and PACF plots of differencing data series were beneficial to estimate the tentative order of the ARIMA model.(4) ARIMA models order obtained based on the evaluation criteria (AIC and BIC) were found skillful.
The residual analysis reinforced the idea of selection of the fitted model order.(5) ARIMA estimates of groundwater storage anomalies fit reasonably well with the observed values as supported by RMSE and NSE skills during the historical training and testing periods.(6) The ARIMA forecasts indicated the increase in March, April and May groundwater storage within the major number of grids of upper CRB.This can be attributed to the early snowmelts in the region during these months as a result of climate change.(7) The study showed a probable decline in future groundwater storage in lower CRB for all months.
Additionally, the model also predicted the decline in near future groundwater storage within upper CRB for the months except for March, April, and May.
The current research will be beneficial to water resource personnel as the newly tested approach showed skillful results.The study can be tested in different regions to verify whether the model is versatile in different varieties of watershed and hydroclimate.The stochastic modeling of GRACE-derived groundwater storage anomalies can be effectively utilized for forecasting groundwater behavior.This will lead in making policies for sustainable groundwater management.Addition to forecasting the ARIMA model coupled with GRACE data can also be utilized to study past variabilities.This stochastic model can also be utilized to determine future water head and trend in CRB which will facilitate water managers to allocate the resources and resolve excessive groundwater withdrawal.Future trend analysis can also be done using the output of the fitted model.A deterministic method can reproduce the spatial scenario but shows poor performance in terms of forecasting calculation.On the other hand, a stochastic method like ARIMA, can't replicate the spatial scenario but can predict the future trend effectively.So, the future study of incorporating ARIMA with another deterministic method can improve model performance.

21 Figure 1 .
Figure 1.Study area showing the Upper and Lower Colorado River Basin.The grids represent the spatial resolution (1° × 1°) of the utilized GRACE data.

Figure 1 .
Figure 1.Study area showing the Upper and Lower Colorado River Basin.The grids represent the spatial resolution (1 • × 1 • ) of the utilized GRACE data.

Figure
Figure 3a,b represent the monthly time series of GRACE-derived ∆TWS, ∆SM, ∆SWE, and ∆WCS, after taking the average from four GLDAS LSMs, from January 2003 to December 2016 for the both Upper and Lower Colorado River Basin.Equation (3) was used to estimate the monthly distributed groundwater storage anomalies based on GRACE using the above-mentioned hydrological components.GRACE-derived ∆GWS was correlated with the average of in-situ measurements in both upper and lower CRB with the R 2 value of 0.62 and 0.67 respectively.The R 2 values testify good conformity between satellite and ground-based measurement as for such analysis the R 2 between 0.55 and 0.75 is ranked as good correlation[21].Like previous studies, the larger study area leads to getting better co-relation among GRACE and in-situ observations.From Figure3a, it can be observed that the groundwater storage in the upper basin declined sharply during 2013.This can be attributed to the lowest recorded snowfall in Rocky Mountain and extreme drought scenario in 2012[52].In lower CRB as observed in Figure3b, groundwater storage followed a decreasing trend from 2004 to 2012 before it experienced a significant depletion from 2012 to mid-2014.On the other hand, the rise in groundwater storage was visible for a shorter period between mid-2009 to mid-2010).Moderately wetted weather lessened the demand for surface water supplies and groundwater recharged during this period.The consistency in groundwater head, in both upper and lower basin, after January 2015 shows strong evidence of recovery from the drought starting from 2000 up to 2014.

Figure 3 .
Figure 3. Average monthly time series of each hydrological component for the period of January 2003-December 2016 for (a) Upper Colorado River Basin and (b) Lower Colorado River Basin; (c) Spatial variation of Thiel Sen Slope at the significance level of p ≤ 0.05.The significance of the trend was obtained with the Mann-Kendal test.

Figure 3 .
Figure 3. Average monthly time series of each hydrological component for the period of January 2003-December 2016 for (a) Upper Colorado River Basin and (b) Lower Colorado River Basin; (c) Spatial variation of Thiel Sen Slope at the significance level of p ≤ 0.05.The significance of the trend was obtained with the Mann-Kendal test.
,b.Two dotted horizontal red lines in the figure signify the margins at 95% CI for calculated autocorrelation and partial autocorrelation values.Here the horizontal axis (x-axis) shows the lagged time steps while the vertical axis denotes the correlation values for the corresponding lags.As seen in the figure, all the correlation values lie between +1 and −1, the maximum and minimum value respectively.Hydrology 2018, 5, x FOR PEER REVIEW 11 of 21 before performing any further analysis.First, the intermediate sample results of ARIMA forecast corresponding to the grid with centroid 110.5 ° W and 31.5 ° N referred hereafter as sample grid is presented.After presenting the results corresponding to the sample grid the results pertaining to all grids of CRB are presented.The ACF and PACF correlograms are drawn in order to determine the stationarity in monthly ΔGWS data as shown in Figure 4a,b.Two dotted horizontal red lines in the figure signify the margins at 95% CI for calculated autocorrelation and partial autocorrelation values.Here the horizontal axis (x-axis) shows the lagged time steps while the vertical axis denotes the correlation values for the corresponding lags.As seen in the figure, all the correlation values lie between +1 and −1, the maximum and minimum value respectively.

Figure 4 .
Figure 4. (a) Sample ACF (b) partial autocorrelation function (PACF) of Gravity Recovery and Climate Experiment (GRACE)-derived ΔGWS (groundwater storage) data for an arbitrary grid and (c) Sample seasonally differenced GRACE-derived ΔGWS data for the sample grid.

Figure 4 .
Figure 4. (a) Sample ACF (b) partial autocorrelation function (PACF) of Gravity Recovery and Climate Experiment (GRACE)-derived ∆GWS (groundwater storage) data for an arbitrary grid and (c) Sample seasonally differenced GRACE-derived ∆GWS data for the sample grid.

Figure 6
Figure 6 depicts the time series plots of ARIMA simulations along with the GRACE-derived monthly ∆GWS for the sample grid.The model training period (January 2003 to December 2014) clearly illustrates the periodicity of existing GRACE-derived ∆GWS.As shown in Figure6, the ARIMA model forecasts were close to the GRACE-derived ∆GWS during the testing period (January 2015-December 2016) signifying good model skills.The model also captured the seasonality and the trend of ∆GWS series remarkably well during the testing period as shown in Figure6.In the training period, the simulated ∆GWS followed the decreasing trend of 0.43 cm/year which was close to the decreasing trend of GRACE-derived ∆GWS evaluated as 0.47 cm/year.On the other hand, during the testing period, both ARIMA results and GRACE-derived ∆GWS showed the recharge in groundwater.Tillman et al.[71] used Coupled Model Intercomparison Project Phase 5 climate projections and supported the idea of recharge in groundwater storage of the CRB.It was observed that the model captured most of the peaks with some minor exceptions.The incompetence of the model to capture some of the peaks may be attributed to the presence of randomness component in the stochastic model.The selected seasonal ARIMA model generates the groundwater storage variability with higher accuracy for the remaining time series.The higher value of Nash-Sutcliffe efficiency (NSE) (0.89) and comparatively lower Root Mean Square Error (RMSE) (1.04) during the training period and 0.72, 1.5 during testing period respectively, indicating higher conformity between satellite-derived and simulated groundwater anomaly for the sample grid.This also suggests robustness of the seasonal ARIMA model that offer anticipated consistency and accuracy.Finally, the tested ARIMA model was utilized to forecast future groundwater anomaly of the sample grid for 36 months.

Figure 6 .
Figure 6.depicts the time series plots of ARIMA simulations along with the GRACE-derived monthly ΔGWS for the sample grid.The model training period (January 2003 to December 2014) clearly illustrates the periodicity of existing GRACE-derived ΔGWS.As shown in Figure6, the ARIMA model forecasts were close to the GRACE-derived ΔGWS during the testing period (January 2015-December 2016) signifying good model skills.The model also captured the seasonality and the trend of ΔGWS series remarkably well during the testing period as shown in Figure6.In the training period, the simulated ΔGWS followed the decreasing trend of 0.43 cm/year which was close to the decreasing trend of GRACE-derived ΔGWS evaluated as 0.47 cm/year.On the other hand, during the testing period, both ARIMA results and GRACE-derived ΔGWS showed the recharge in groundwater.Tillman et al.[71] used Coupled Model Intercomparison Project Phase 5 climate projections and supported the idea of recharge in groundwater storage of the CRB.It was observed that the model captured most of the peaks with some minor exceptions.The incompetence of the model to capture some of the peaks may be attributed to the presence of randomness component in the stochastic model.The selected seasonal ARIMA model generates the groundwater storage variability with higher accuracy for the remaining time series.The higher value of Nash-Sutcliffe efficiency (NSE) (0.89) and comparatively lower Root Mean Square Error (RMSE) (1.04) during the training period and 0.72, 1.5 during testing period respectively, indicating higher conformity between satellite-derived and simulated groundwater anomaly for the sample grid.This also suggests robustness of the seasonal ARIMA model that offer anticipated consistency and accuracy.Finally, the tested ARIMA model was utilized to forecast future groundwater anomaly of the sample grid for 36 months.

Figure 6 .
Figure 6.GRACE-derived ΔGWS and ARIMA resulted with the best fit model values and forecasted ΔGWS.

Figure 6 .
Figure 6.GRACE-derived ∆GWS and ARIMA resulted with the best fit model values and forecasted ∆GWS.

Figure 7 .
Figure 7. ARIMA model order for each grid cells representing (a) autoregressive term (b) differencing term and (c) moving average term for both (i) non-seasonal and (ii) seasonal component for the Colorado River Basin.
summarizes the distribution of RMSE and NSE respectively during the training and testing period for CRB.The RMSE and NSE were obtained from the simulated and GRACE-derived ΔGWS data presented in Figure 8a,b.The red line within the box plots of Figure 8c,d represent the median value while the horizontal box edges represent the 25th and 75th percentile.The whiskers represent the 5th and 95th percentile.Further,

Figure 7 .
Figure 7. ARIMA model order for each grid cells representing (a) autoregressive term (b) differencing term and (c) moving average term for both (i) non-seasonal and (ii) seasonal component for the Colorado River Basin.
summarizes the distribution of RMSE and NSE respectively during the training and testing period for CRB.The RMSE and NSE were obtained from the simulated and GRACE-derived ∆GWS data presented in Figure 8a,b.The red line within the box plots of Figure 8c,d represent the median value while the horizontal box edges represent the 25th and 75th percentile.The whiskers represent the 5th and 95th percentile.Further, the robustness of the ARIMA model is established by NSE during both training and testing period which varies between 0.5 and 1.

Figure 9 .
Figure 9. Box plot of the forecasted monthly groundwater anomalies of (a) Upper Colorado River basin and (b) Lower Colorado River basin.

Figure 9 .
Figure 9. Box plot of the forecasted monthly groundwater anomalies of (a) Upper Colorado River basin and (b) Lower Colorado River basin.
Hydrology 2018, 5, x FOR PEER REVIEW 12 of 21 close to the tentative order to estimate ARIMA model parameters.The comparison of AIC and BIC of different models are summarized in Table

Table 1 .
AIC and BIC values for different models.

Table 1 .
AIC and BIC values for different models.