GURU PRADHAN March, 2014
ITC SUPERVISOR IIRS SUPERVISORS Dr. Rogier van der Velde Dr. P.K. Champati Ray Mr. Suresh Kannaujiya
Understanding
Interannual Groundwater
Variability in North India
using GRACE
Thesis submitted to the Faculty of Geo-information Science and Earth Observation of the University of Twente in partial fulfilment of the requirements for the degree of Master of Science in Geo-information Science and Earth Observation.
Specialization: Natural Hazards and Disaster Risk Management
THESIS ASSESSMENT BOARD:
Chairperson : Dr. V.G. Jetten
External Examiner : Dr. Rajesh S., (WIHG, Dehradun) ITC Supervisor : Dr. Rogier van der Velde
IIRS Supervisor : Dr. P.K. Champati Ray IIRS Supervisor : Mr. Suresh Kannaujiya
OBSERVERS:
ITC Observer : Dr. N.A.S. Hamm
Interannual Groundwater Variability in North India using GRACE
Guru Pradhan
Enschede, the Netherlands [March, 2014]
DISCLAIMER
This document describes work undertaken as part of a programme of study at the Faculty of Geo-information Science and Earth Observation (ITC), University of Twente, The Netherlands. All views and opinions expressed therein remain the sole responsibility of the author, and do not necessarily represent those of the institute.
“When the well is dry, we know the worth of water.”
- Benjamin Franklin
“For all of its uncertainty, we cannot flee the future”
- Barbara Jordan
“Uncertainty and mystery are energies of life. Don’t let them scare you unduly,
for they keep boredom at bay and spark creativity.”
- R.I. Fitzhenry
ABSTRACT
Though groundwater stress in North India has been extensively studied in the past with the help of the Gravity Recovery and Climate Experiment (GRACE) satellite mission, interannual variability of groundwater storage (GWS) in the region remains unexplored. Higher variability denotes more uncertainty and volatility, and an uncertain groundwater supply can lead to unexpected shortfalls during periods when water is needed the most. From January 2003 to December 2012, GWS depletion rate was found to be approximately -1.6 ± 0.04 cm yr-1 in the area roughly corresponding to the Indian side of the Indo- Gangetic plains which translates into about 160 km3 of net groundwater loss – equivalent to more than sixteen times the active capacity of India’s largest dam reservoir. In the Indus and Ganges Basins, the groundwater depletion rate has slowed down remarkably compared to previous studies but there are still respective net groundwater losses of 60 km3 and 100 km3 over the 10-year study period. Furthermore, groundwater depletion was found to be particularly pronounced near the Himalayan inter-plate collision zone suggesting leakage of tectonic and erosion-driven mass loss signals into the GWS solution.
The interannual standard deviation of groundwater was found to be 1.5 ± 0.1 cm, which translates into a projected 95% confident yearly regional groundwater shortfall of 34 km3. This unexpected deficit can be particularly devastating during a dry year. Moreover, the Gangetic basin was found to have a higher interannual standard deviation than the Indus Basin indicating higher groundwater variability in the former. Subsequently, this work went on to study the time evolution of interannual GWS variability and detected a non-negligible rise in the 12-month moving standard deviation over North India of 0.12 ± 0.04 cm yr-1. This suggests that regional GWS variability is on the rise. The interpretation of this development is difficult but increased monitoring of groundwater in the region is warranted as this change might be a precursor to a fundamental shift in the regional groundwater system.
The GRACE-derived GWS values agree reasonably well with in situ well observations with correlation of 0.80 ± 0.03 and R-Squared of 0.64 ± 0.05 over a part of Uttar Pradesh thus supporting the validity of this remote sensing approach. The RMS error between the GRACE and in situ GWS time series was 8.43 ± 0.53 cm which is quite high however considering inherent uncertainties in the in situ data this study considers this margin of error reasonable. This work concludes that rapid, accurate regional GWS mapping is possible through the GRACE mission. In light of continued groundwater depletion and possible increase in GWS variability over North India, this study advocates for more robust water conservation and storage measures.
Keywords: GRACE, GLDAS, Groundwater, Northern India, Interannual Variability, CGWB, Depletion, Water Security, Hydrology, Well, Empirical Bayesian Kriging
ACKNOWLEDGEMENTS
This work could not have been possible without the help of many sources – both human and digital!
Countless people have contributed directly or indirectly in the completion of this thesis and I thank you all from the bottom of my heart! Please do not feel offended if you don’t find your name explicitly mentioned here – you are still important!
I would like to express my deepest gratitude to Dr. Y.V.N. Krishna Murthy (Director, IIRS), Dr. S.P.S.
Kushwaha (Group Director, PPEG), Dr. P.L.N Raju (Group Head, Geoinformatics), Dr. S.K. Saha (Group Head, Earth Resources and Systems) for facilitating my transition from Post-Graduate Diploma to the Masters of Science program. Countless thanks to my supervisors in Indian Institute of Remote Sensing and University of Twente – Dr. P.K. Champati Ray, Dr. Rogier van der Velde, Mr. Suresh Kannaujiya – for the constant encouragement and crisp, indispensable feedback.
Of course, I am indebted to the faculty and research fellows from each academic department at IIRS who gave me a brief but insightful introduction to topics and matters in fields ranging from marine &
atmospheric systems to urban information systems. I also extend my sincere gratitude to the instructors and research fellows who enlightened me in the advanced modules and coursework at University of Twente.
Now, a word of appreciation to my classmates, seniors, and friends at IIRS and UT; there are far too many of you to mention but I will try anyhow. Raj B., Colonel Sir, Ishaan, Kanishk, Aravind, Piyush, Abhisheikh S., Tanya, Akhil, Apoorva, Neha, Hemlata, Ravisha, Amit, Shreya, Girija, Sahil, Pradeep – it is a good MSc. group we got, we will stay in touch. Unmesh, Parag, Amit Sir, Shashi G., Anant, Hemanthi, Ponraj, Sai K., Bond, Shishant, Soumya, Joyson, Akanksha, Vivek – my MTech. comrades, it was a pleasure. Special mention goes to Rohit, the best roommate out there, and Taibang for his gracious hospitality and expert culinary skills. My long-departed PG Diploma friends – Kelhou, Priyanka, Rosmary, the Rafaels, Oswaldo, David, Ramanuj, Vedika, Ritanjali, Manish Sir, Gangesh, Himanshu, Tarun, Moa, Azzedine, Tamrat, Preeti, Anil, Blecy – miss you all and hope you are doing well. My ITEC compatriots – Dorji, Daniel, Jasneel, Menuka, Angela, Sabin, Calvin, Sherzhod, Thabo, Eduardo, Ganessan, Ajit, Tia – we shall meet again. Friends at UT – Zelalem, Raghudeep, Iria, Twang, Franklin, Payam, Mehdi, Sidhi, Gunamani, Rahul Raj – it was fun!
Space is short so I will list the remaining folks here. Thank you all! My parents and extended family, librarians, administrative staff, mess members, security staff, cleaning staff, Lucy the campus pet dog.
TABLE OF CONTENTS
List of Figures... V List of Tables ... VI
1. INTRODUCTION ... 1
1.1. Background ... 1
1.2. Related Work ... 2
1.3. Research Identification... 4
1.3.1. General Objective ... 4
1.3.2. Sub-Objectives ... 4
1.3.3. Research Questions ... 4
1.4. Thesis Structure ... 5
2. THEORETICAL BACKGROUND ... 7
2.1. Remote Sensing Tools and Land Surface Modelling ... 7
2.1.1. GRACE Mission ... 7
2.1.2. GLDAS ... 9
2.2. In Situ Well Data ... 9
2.3. Estimation of Groundwater Storage Change Using Water Balance Approach ... 10
2.4. Interannual Variability of Groundwater Storage ... 11
2.4.1. Interannual Standard Deviation ... 11
2.4.2. Shift in Interannual Groundwater Storage Variability ... 12
2.5. Error Propagation ... 13
2.5.1. First-Order Error Propagation ... 14
2.5.2. Monte Carlo Method ... 15
3. STUDY AREA ... 16
3.1. Background ... 16
3.2. Agriculture ... 17
3.3. Climate ... 18
3.4. Hydrogeology ... 18
4. DATA AND METHODS ... 20
4.1. GRACE-Derived Terrestrial Water Storage ... 20
4.1.1. Land Grid Scaling ... 20
4.1.2. Error Handling ... 20
4.2. GLDAS Land Surface Models ... 21
4.2.1. Anomaly and Uncertainty Expression ... 21
4.3. In Situ Well Data ... 23
4.3.1. Quality Control Mechanism ... 23
4.3.2. Well Time Series Resampling ... 25
4.3.3. Selection of Validation Region ... 25
4.3.4. Computation of Groundwater Storage Fluctuations ... 25
4.3.5. Gridding of Groundwater Storage Anomalies Using Empirical Bayesian Kriging ... 27
4.3.6. Cross-Validation and Prediction Errors ... 30
4.4. Software Used ... 30
4.5. Methodology ... 30
4.5.1. Study Area Masking... 31
4.5.2. Deseasonalization ... 32
4.5.3. Trend Analysis ... 33
4.5.4. Detrending ... 33
4.5.5. Regional Averaging ... 33
4.5.6. In Situ Validation ... 34
4.5.7. Error Analysis ... 34
5. RESULTS AND DISCUSSION ... 36
5.1. Groundwater Storage Time Series ... 36
5.1.1. Indo-Gangetic Plains ... 36
5.1.2. Northwest India ... 39
5.1.3. Ganga Basin ... 40
5.1.4. Discussion ... 41
5.2. Interannual Standard Deviation ... 41
5.2.1. Indo-Gangetic Plains ... 41
5.2.2. Northwest India ... 42
5.2.3. Ganga Basin ... 42
5.2.4. Discussion ... 42
5.3. Shifts in Interannual GWS Variability ... 43
5.3.1. Indo-Gangetic Plains ... 43
5.3.2. Northwest India ... 44
5.3.3. Ganga Basin ... 44
5.3.4. Discussion ... 44
5.4. Trend Analysis Comparison and Variability Statistics ... 45
5.5. Selected Validation of GRACE-Derived GWS in Uttar Pradesh Sub-Region ... 46
6. CONCLUSION AND RECOMMENDATIONS ... 49
6.1. Conclusion ... 49
6.2. Recommendations ... 52
REFERENCES ... 54
APPENDICES... 62
Appendix - I: Land Surface Model Parameters ... 62
Appendix - II: Jarque-Bera Test ... 64
Appendix - III: Empirical Bayesian Kriging Parameters ... 65
Appendix - IV: Kriging Cross-Validation Results ... 66
LIST OF FIGURES
Figure 2-1: The GRACE Mission ... 8
Figure 2-2: Visualization of equivalent water thickness (EWT) ... 9
Figure 2-3: Schematic for visualizing effects on extreme GWS events in different scenarios .. 12
Figure 3-1: Political map of India with the study area states labelled... 16
Figure 3-2: Agriculture zone map of Indo-Gangetic Plains . ... 17
Figure 3-3: Hydrogeological Map of India ... 18
Figure 4-1: Well Data Processing Workflow... 24
Figure 4-2: Location of Short-Listed Wells ... 24
Figure 4-3: Flow chart describing overall methodology ... 31
Figure 4-4: Study Area Mask ... 32
Figure 5-1: Indo-Gangetic Plains GWS Time Series ... 36
Figure 5-2: Deseasonalized Trend Map of Indo-Gangetic Region ... 38
Figure 5-3: Rodell and Tiwari Study Mask ... 38
Figure 5-4: Indus Basin Study Mask, Raw & Deseasonalized GWS Plots. ... 39
Figure 5-5: Ganges Basin Study Mask, Raw & Deseasonalized GWS Plots. ... 40
Figure 5-6: Interannual Standard Deviation Map of Indo-Gangetic Region ... 41
Figure 5-7: Moving Statistics for Indo-Gangetic Plains ... 43
Figure 5-8: Moving Statistics for Northwest India ... 44
Figure 5-9: Moving Statistics for Ganga Basin ... 44
Figure 5-10: Validation over Uttar Pradesh ... 46
LIST OF TABLES
Table 2-1: First Order Error Propagation Approximation Formulae ... 14
Table 4-1: Primary and Secondary Source Datasets ... 22
Table 4-2: List of Uncertainties When Calculating GWS for In Situ Well Data ... 27
Table 5-1: Trend Analysis Comparison between Present and Previous Works ... 45
Table 5-2: Regional Statistics Table ... 46
Table 5-3: Validation Statistics ... 47
Table A: Physical parameters available from GLDAS ... 61
Table B: Soil parametrization of the land surface models ... 62
Table C: Error Statistics of Kriging Process ... 65
1. INTRODUCTION
1.1. Background
North India has earned the unenviable distinction of being one of the most water stressed regions in the world (“India among high risk nations in water stress survey,” 2013). With a population exceeding 500 million people, the region’s agriculture and food security is heavily reliant on groundwater-based irrigation to feed its large yet booming population. The combination of climate change, large-scale contamination of shallow groundwater resources (Nahar, Hossain, & Hossain, 2008), and the planned diversion of surface water resources from the Northern regions to the dryer Southern regions (Misra et al., 2006) all serve to intensify an already deteriorating water security situation. Furthermore, recent regional groundwater mapping of North India using Gravity Recovery and Climatology Experiment (GRACE) satellite observations and Global Data Assimilation System land surface models (see Chapter 2 for more detailed discussion regarding measurement of GWS) have all unanimously shown unsustainable groundwater depletion trends (Gleeson et al., 2012; Rodell et al., 2009; Tiwari et al., 2009). Though these findings helped quantify and illuminate the worsening groundwater scenario, they have not fully exploited the immense information contained within the GRACE-derived groundwater storage (GWS) time series.
This study has multiple overlapping goals which ultimately link to a better understanding of the interannual groundwater variability over North India from 2003 till 2012. The first goal is to re-calculate and update the groundwater depletion trend (if any) across the North Indian region to assess any shifts or changes in regional groundwater activity with respect to GRACE estimates from previous studies. The work then moves beyond the GWS trend analysis and explores the interannual variability of groundwater over North India. This will be achieved by first quantifying the level of year-to-year variability across the region and then studying the time evolution of yearly GWS variability to identify possible shifts in regional groundwater dynamics over the 10-year study period. Another novel aspect of this research is the validation of in situ well data with GRACE observations over a part of the Indian state of Uttar Pradesh.
Interannual variation of GWS measures the variability of groundwater levels at time scales of one (1) year or longer. Groundwater fluctuations are driven by natural (rainfall, vegetation, soil types) and anthropogenic (socio-economic concerns, land use/land cover change, damming) processes with complex, non-linear interactions between them. Thus, one can expect to see a wide range of variability in yearly groundwater storage with some areas facing unexpected shortages or flooding. As such, interannual variability of GWS is a gauge of groundwater supply unpredictability and should be combined with mean GWS depletion rate to comprehensively measure water stress. Areas experiencing high interannual GWS variability face higher risk of water supply shortages even if there is little to no net loss of groundwater. In
the absence of adequate water governance and storage mechanisms, these supply shocks can be especially devastating during dry years (Reig et al., 2013). The level of yearly GWS variability can be captured by simply calculating the standard deviation of yearly average GWS in a manner similar to the one used to calculate interannual variation of sunlight over United States (Wilcox & Gueymard, 2010). More details on this can be found in Chapter 2.2.
This study also tackles the time evolution of GWS interannual variability over North India. Inspired by studies of shifts in river discharge variability, this work uses the moving standard deviation as the primary method to analyse the progression of yearly GWS variability (Booth et al., 2006). Even in areas where there is no net GWS change, any increase in GWS variability translates into higher frequency of extreme GWS behaviour (both high and low water tables), as well as an increase in the magnitude of these extreme events. Furthermore, a shift in the variability may be a precursor to fundamental changes in the groundwater – and perhaps underlying large-scale hydrological – dynamics which could lead to unexpected challenges that the region must contend with. An increase in GWS variability, combined with a negative trend in GWS levels, is even more worrisome. Such a development would denote lower water tables with more extreme drops in the water level which can have harsh implications for groundwater access and supply. The GWS depletion rates calculated during the course of this study will be combined with interannual variability results to better assess regional groundwater behaviour.
The final problem this thesis tackles involves the validation of the GRACE-derived GWS solution with well observations provided by Central Groundwater Board (CGWB), the primary authority for groundwater surveying and monitoring in India. Indeed, robust validation of the GRACE solution further reinforces the product’s legitimacy for regional groundwater mapping applications. Due to both temporal and spatial undersampling of well observations, this work only considers selective validation of in situ data over the state of Uttar Pradesh. This research uses Empirical Bayesian Kriging (EBK), a relatively new geostatistical interpolation method to interpolate CGWB-derived GWS values. This novel approach corrects for and automates some model fitting difficulties inherent in classical kriging (Krivoruchko, 2012). Using datasets from the GRACE mission, the Global Land Data Assimilation System (GLDAS) land surface models, and CGWB; this thesis is a preliminary study of interannual variability of groundwater storage in the North Indian region.
1.2. Related Work
The regional analysis of groundwater depletion over North India using GRACE has been explored independently by both Tiwari et al. (2009) and Rodell et al. (2009) with both parties concluding that large- scale, unsustainable groundwater extraction is taking place (Rodell et al., 2009; Tiwari et al., 2009). This work was followed by a groundwater stress map or footprint of the Upper Ganges aquifer by Gleeson et
a hydrological model and country/administrative unit water use statistics (Gleeson et al., 2012). All these works do well in characterizing first-order effects such as net depletion rates but they do not incorporate second-order effects such as the variability of the groundwater storage. Variability of groundwater storage translates into supply unpredictability which must be accounted for in water security assessments and met with adequate water storage and governance mechanisms (Reig et al., 2013).
A qualitative analysis of interannual variability of terrestrial water storage using GRACE has been explored over South Asia by tracking changes in the yearly average TWS values (Shum et al., 2010). The aforementioned study, however, does not focus on GWS and the study area encompasses a far larger area.
Interannual variability of groundwater observation well data has been studied in Canada at three different regions using wavelet transforms and has shown remarkably different variability behaviours for each of the sites (Tremblay et al, 2011).
A coefficient of variation map has been prepared in many domains to quantify interannual variability for different phenomena: solar radiation (Wilcox & Gueymard, 2010), river discharge (Booth et al., 2006;
Restrepo & Kjerfve, 2000), and NDVI (Milich & Weiss, 2000). After an extensive literature survey, no maps for interannual variability of GWS were found. The closest thing uncovered was the interannual variability map of total blue water detailed in the Aqueduct Water Risk Map (Reig et al., 2013). However, this map measures the interannual variability of available blue water which is essentially a measure of the runoff flowing into a catchment as measured from GLDAS simulations.
The main method for understanding and quantifying shifts in interannual GWS variability has been borrowed from literature rooted in climate change studies (Folland et al., 2002; Vinnikov & Robock, 2002). Folland et al. (2002) used extensive climatological datasets spanning nearly 100 years of observations and studied the probability distribution of both precipitation and temperature across different time periods to uncover any shifts in climate variability. Vinnikov and Robock, on the other hand, fit trends through the moments of various climatic indices to determine whether the observed climate is getting more or less variable. Vinnikov and Robock’s approach is adopted in this thesis to test for shifts in GWS variability.
Validation of the GRACE-derived groundwater storage solutions have been carried out successfully over areas as diverse as Mississippi River Basin (Rodell et al., 2006), Bangladesh (Shamsudduha et al., 2012), and Yemen (Moore & Fisher, 2012) over study areas exceeding 200,000 km2. As of now, no extensive validation of GRACE has been carried out over India using in situ well observations. Additionally, whenever efforts to validate GRACE results have been carried out, they have all used simple averaging of in situ well data or deterministic interpolation techniques such as inverse distance weighing (IDW) or Thiessen polygons to derive gridded GWS outputs. This is contrary to the rich and diverse literature all
concurring with the effectiveness of geostatistical methods (kriging) for hydrological applications – especially with interpolation of groundwater data (Ahmadi & Sedghamiz, 2006; Delbari et al., 2013;
Machiwal et al., 2012). Even then, these studies all use ordinary kriging methods where several theoretical semivariograms (spherical, exponential, circular etc.) are used to model the actual or empirical semivariogram. The semivariogram model that best fits the empirical semivariogram is then assumed to be the true semivariogram without taking into account fitting errors. This study departs from previous kriging efforts by using Empirical Bayesian Kriging in order to correct for this fitting problem.
1.3. Research Identification
1.3.1. General Objective
This study aims at quantifying and understanding interannual variability of groundwater storage over the North Indian states of Bihar, Haryana (including Delhi NCR), Punjab, Rajasthan, Uttar Pradesh, and West Bengal for the time period encompassing January 2003 till December 2012. Better understanding of interannual GWS variability would require derivation of groundwater depletion rates (if any) and also validation of GRACE data with in situ well observations to bolster the legitimacy of this remote sensing approach.
1.3.2. Sub-Objectives
The successful accomplishment of the research objective requires that the following sub-objectives be met:
To estimate the GWS and its rate of change in the North Indian region from 2003 till 2012, and compare with previous works on the same subject
To compute the interannual standard deviation both at the pixel (1° x 1°) and regional level in order to quantify interannual GWS variability
To detect any change in interannual GWS variability in North India during the period 2003- 2012
To validate in situ well observation data with GRACE-derived GWS solution for the state of Uttar Pradesh using Empirical Bayesian Kriging
1.3.3. Research Questions
The study will attempt to answer the following research questions:
1) How is GWS estimated from GRACE and GLDAS observations?
2) What is the level of uncertainty associated with the GRACE-derived GWS results and how do these errors propagate?
3) How is GWS estimated from in situ well observations?
4) What criteria are to be used to enable effective comparison between GRACE and CGWB- derived GWS solutions?
5) What can the interannual standard deviation say about groundwater dynamics in the study region?
6) Can information regarding GWS depletion rate and time evolution of groundwater variability be used to shed light on regional water stress, and what might be effective ways to combat it?
1.4. Thesis Structure
The research work is organized as follows:
Chapter 1: Introduction – the concept of groundwater interannual variability is introduced and the thesis’s motivation is described. The research objectives of this study and the research questions that are to be answered are furthermore presented along with previous work carried out.
Chapter 2: Theoretical Background – this section deals with the derivation of GWS using the Water Balance Method from GRACE Terrestrial Water Storage (TWS) measurements and GLDAS hydrological outputs.
The concept of interannual variability and measurement of shifts in variability is then displayed.
Subsequently, there is an introduction to the GRACE satellite mission, the GLDAS project, and the CGWB well data. Then, the importance of error and uncertainty propagation is expanded upon, and the use of Monte Carlo methods is commented upon. The chapter ends with the literature review.
Chapter 3: Study Area – the climatic, hydrogeological, agro-economic context of the North Indian region is dealt with in this chapter.
Chapter 4: Data and Methods – this chapter deals with the datasets used in the course of this research. A short description regarding the GRACE, GLDAS, and CGWB datasets is given and is followed by the software used. Various techniques utilized for processing GRACE, GLDAS, CGWB datasets are explained in detail: Trend Analysis, Time Series Decomposition, Monte Carlo Method, Calculation of Interannual Standard Deviation, Moving Statistics, Regional Analysis, Well Data Processing, and Empirical Bayesian Kriging.
Chapter 5: Results and Discussion - the interannual variation level across North India is discussed and then followed by results regarding changes in interannual GWS variability. Subsequently, the validation efforts are analysed and the accuracy of GRACE is commented upon.
Chapter 6: Conclusion and Recommendations – the final section provides conclusions and recommendations for future work along with some advice for policymakers.
2. THEORETICAL BACKGROUND
2.1. Remote Sensing Tools and Land Surface Modelling
Answering the research questions required the use of the GRACE Level 3 Processed Terrestrial Water Storage (TWS) solutions and the land surface model outputs from GLDAS. A short description of both these datasets is outlined below.
2.1.1. GRACE Mission
The Gravity Recovery and Climate Experiment (GRACE) is unprecedented in that it is the first satellite remote sensing mission directly applicable for regional groundwater mapping (Rodell et al., 2006) though its primary directive is to obtain accurate estimates of Earth’s gravity field variations (Tapley et al., 2004). The mission consists of two satellites flying in tandem in a polar, near circular orbit at 500 km altitude with an inter-satellite separation distance of approximately 220 km. The gravity field information is actually inferred from the inter-satellite distance which is measured within µm accuracy using a K-Band microwave system (Tapley et al., 2004). Potential error sources such as atmospheric drag and satellite perturbations are measured and filtered out using readings from a highly accurate, on-board accelerometer while precise positioning is determined using on-board GPS receivers (Tapley et al., 2004).
The principal and admittedly unconventional idea behind GRACE is that the satellites themselves act as the principal measurement devices. The two satellites (also known as ‘Tom’ and ‘Jerry’) work in tandem to map the gravitational field of the Earth. When surface features that distort the gravitational strength such as mountains (which decrease the gravitational field strength) are encountered, the leading satellite accelerates by a certain amount followed by the trailing satellite which then catches up (see figure 2.1). These minute changes in the inter-satellite distance are then fed into what are essentially massive regression engines in order to determine the gravity field strength at the data processing facilities in Jet Propulsion Laboratory (JPL), University of Texas – Austin Centre for Space Research (CSR), and the GeoForschungZentrum (GFZ).
Figure 2-1: The inter-satellite distance varies according to the local gravitational anomaly. When the leading satellite encounters a land/sea interface (left), it decelerates due to the higher gravitational pull because of the denser oceanic crust. On (right), mountainous regions have lower gravitational strength
due to isostasy so the inter-satellite distance increases as the leading satellite accelerates (Haile, 2011).
At the data centres, the K-Band Level-1 radar data is converted into gravitational field information in the form of spherical harmonic coefficients (Rodell et al., 2006; Wahr et al., 1998). These coefficients are available as Level-2 products and have been corrected for atmospheric and oceanic circulation, and solid Earth tides using underlying models (Rodell et al., 2006).
Water is quite heavy and is usually the largest contributor to mass variability on the earth’s surface (Ogawa, 2010). By assuming that gravity changes are primarily driven by changes in distribution of terrestrial water storage, the Level-2 spherical harmonics data can be further processed to express the gravity field changes in terms of ‘equivalent water thickness’ (Ogawa, 2010). The term ‘equivalent water thickness’, or EWT, is derived from the assumption that these hydrological mass changes are concentrated in a very thin layer on the Earth’s surface (Wahr et al., 1998) whose vertical extent is measured in centimetres (see fig 2.2). Changes in EWT are available as Level-3 Terrestrial Water Storage (TWS) solutions after further filtering and correction for post-glacial rebound. This dataset provides the values needed for the water balance equation detailed in equation (2.1) and (2.2). Furthermore, this processing limits the accuracy of TWS observations to be valid only for large geographical areas exceeding 200,000 km2.
Figure 2-2: Visualization of equivalent water thickness (EWT). Note that the terrestrial water storage is assumed to be concentrated in a very thin layer on the Earth’s surface (Ogawa, 2010)
2.1.2. GLDAS
The Global Land Data Assimilation System (GLDAS) is an inter-institutional effort undertaken by National Aeronautics and Space Administration (NASA) and National Oceanic and Atmospheric Administration (NOAA) in order to simulate land surface state (soil moisture, surface temperature) and flux (evapotranspiration, sensible heat flux) parameters globally at high spatial resolution and in a near real-time basis (Rodell et al., 2004). As such, it drives four (4) offline (uncoupled with atmosphere) land surface models, assimilates an enormous amount of satellite and ground-based observations, and produces outputs at resolutions ranging from 0.25° to 1°. The land surface models being driven by GLDAS are: NOAH, MOSAIC, Community Land Model (CLM), and Variable Infiltration Capacity (VIC). The principal advantage of GLDAS is its sophisticated data assimilation engine which ingests vast amounts of remote sensing and in situ observations in a near real-time fashion in order to constrain land surface model outputs from deviating too far from observed states. The components of the water balance equation detailed in section (2.3) of this chapter are obtained using GLDAS land surface model simulations.
2.2. In Situ Well Data
This study has access to data from 3597 monitoring wells scattered across the Indian states of Delhi NCR, Haryana, Rajasthan, and Uttar Pradesh kindly provided to us by Central Groundwater Board (CGWB). Initially, most of the CGWB well observation system consisted of dug wells. However, these are being gradually replaced by piezometers for water level monitoring with measurements generally
being taken four times a year in the months of January, April/May, August, and November (Jain &
Singh, 2003). Recently, automated measurement systems are being installed in select wells that would enable real-time monitoring of groundwater levels however nationwide implementation is yet to take place (Goswami, 2014).
Even though the guidelines stipulated quarterly observation of groundwater levels, the periodicity of these measurements vary state to state, and from site to site (Central Ground Water Board, 2011).
Accordingly, the CGWB datasets were subject to systematic quality control procedures to screen out well records with large temporal gaps, missing records, or insufficient observations. Subsequently, the short-listed well time series were converted into groundwater storage changes, and then gridded to the same spatial resolution as the GRACE Level-3 TWS solution. Furthermore, the CGWB well data was processed only for the state of Uttar Pradesh in order to selectively validate the GRACE solutions. A detailed description and rationale of this methodology is described in Chapter 4.
2.3. Estimation of Groundwater Storage Change Using Water Balance Approach
Estimation of groundwater storage (GWS) at a regional scale can be approximated through the use of the water balance equation. A water balance approach is a fundamental hydrological technique which states that the flow of water into and outside of a system must equalize or ‘balance’ (Rodell et al., 2006) with the change in storage. Accordingly, the water balance method to compute terrestrial water storage changes in the North Indian region can be expressed as:
(2.1)
Where is the change in terrestrial water storage, is the change in soil moisture, is the change in groundwater storage, is the change in snow-water equivalent, represents the change in canopy storage, and is the change in surface water storage. All the above parameters are expressed as centimetres of equivalent water thickness (EWT). Re-arranging equation (2.1) in context of the GRACE-derived Level 3 (L3) estimates, the GLDAS-based land surface model outputs of soil moisture , snow-water equivalent , and canopy storage , and assuming that surface water change is negligible; we can calculate change in groundwater storage as:
(2.2)
The GRACE L3 solution vertically integrates all sources of hydrological variability as well as other unaccounted sources such as earthquakes (Mikhailov et al., 2004) and other residuals. The decision to
discount surface water change stems from the general observation that surface water can be considered as the intersection of the water table with the land surface (Winter, 1999). Since the approach taken in equation (2.2) assumes that groundwater is spatially continuous across the area of interest, surface water can be considered as an extension of groundwater and can be removed from the water balance equation.
Still, this assumption fails during times of extreme flooding and in very moist regions of the world such as the Amazon (Rodell et al., 2006). Nevertheless, as a general approximation for GWS, equation (2.2) holds value and is used in this thesis.
2.4. Interannual Variability of Groundwater Storage
Studying GWS variability is often confounded by the presence of seasonal and trend components. These factors often exaggerate or dampen the actual variability of the time series so they must be removed in order to properly isolate the variability of the groundwater levels (Zhang & Qi, 2005). The raw GWS time series must be de-seasonalized and detrended (see section 4.5) in order to isolate the groundwater storage variations. Subsequently, the task of quantifying and detecting shifts in GWS interannual variability can be undertaken.
2.4.1. Interannual Standard Deviation
In order to quantify the level of year-to-year volatility in GWS levels, the interannual standard deviation is computed both at the pixel and the regional level. The interannual standard deviation is derived in a manner similar to Wilcox & Gueymard, (2010) and will retrieve the mean annual GWS over the 10-year period and the average annual GWS value for each specific year to calculate the interannual standard deviation using the following formula:
√[ ∑ ( ) ] (2.3)
The result of equation (2.3) will express the general level of dispersion around the mean yearly, long- term groundwater storage level. It is a measure of the volatility of the groundwater supply and conveys the likelihood of data falling within around the mean. Assuming a normal distribution of GWS, one can be reasonably or 68% confident that the mean yearly GWS level will be within one around the long-term mean. However, when designing robust water management policies it makes sense to work with 95% confidence intervals and make allowances that yearly mean GWS levels may fluctuate by from the long-term average. The higher the interannual standard deviation, the larger the likelihood of extreme swings in GWS levels on a year-to-year basis. This is, of course, a simplification as groundwater varies in complex ways in response to external stimuli however the interannual
standard deviation can serve as a useful approximation of yearly GWS volatility and as a rough tool for making yearly groundwater shortfall projections.
2.4.2. Shift in Interannual Groundwater Storage Variability
It is important to differentiate between changes in the mean GWS levels (which are approximated through trend analysis) and changes in GWS variability. Changes in GWS variability can be thought as the change in the range of the highest and lowest GWS levels (Hiraishi et al., 2000). Thus, an increase in GWS variability translates into an increase in the frequency of both extreme high and low water table fluctuations, as well as an amplification of the magnitude of these events. Though there has been a lot of emphasis on declining groundwater levels, extremely high water tables are just as disruptive. If the water table rises to the surface, one can expect flooding of urban areas, waterlogging and salinization of agricultural fields, and more amenable conditions for growth of pests such as mosquitoes, ticks, locusts, rodents etc. (Kovalevsky, 1992). Thus, it is of immense value that any shifts in GWS variations are understood and accounted for when implementing water storage and distribution schemes. Figure 2-3 demonstrates different scenarios involving changes in the groundwater regime and their effects on the water table.
Figure 2-3: Schematic for visualizing effects on extreme GWS events in different scenarios
Probability of Occurrence
The methodology adopted to identify shifts in GWS interannual variability is quite a simple one and involves the use of moving statistics. Inspired by the use of moving standard deviation for detecting shifts in streamflow variability in the hydrological community (Booth et al., 2006), a simple 12-month moving standard deviation is passed over the GWS time series to identify any shifts in the yearly standard deviation. The moving standard deviation is calculated for an arbitrary window size in the following manner:
̂
{
√ ( ( ̂ ) ( ̂ ) ( ̂ ) )
√ (( ̂ ) ( ̂ ) ( ̂ ) )
(2.4)
Where ̂ is the simple moving average at time and is the last time index value recorded. Similar methods can be extended to calculate the moving skewness and moving kurtosis of a time series. A 12- month running skewness window will track changes in skewness of the dataset and can be used to detect changes in direction of extreme events (Doane & Seward, 2011). For example, an increase in skewness for GWS anomalies suggests a longer right-tail. Incidentally, a right-tailed distribution indicates a larger tendency towards positive GWS values and an uptick in the frequency of higher than normal water tables thus leading to more urban flooding events, waterlogging etc. Similarly, a decreasing skewness suggests a propensity towards negative GWS events as well as in the number of extreme negative GWS anomalies.
A 12-month moving kurtosis window will detect shifts in yearly kurtosis – a statistical measure of how heavy tailed a random variable is (DeCarlo, 1997; “SAS Elementary Statistics Procedure,” 2008). High kurtosis signifies that a large portion of the overall variance is due to a few extreme events. A lower kurtosis, on the other hand, suggests a more uniform distribution with the variance primarily driven by a large number of small- to medium-sized deviations from the mean.
2.5. Error Propagation
No scientific analysis, model, or simulation is complete without accounting for uncertainty in the raw datasets, and its propagation through intermediate and final calculations. The GRACE Level-3 TWS datasets have significant measurement and leakage (from oceans and surrounding pixels) errors associated with each pixel. Different land surface model outputs for the same variable and for the same location may vary considerably due to different modelling approaches, algorithm implementations, and soil layer parameterisations (Rodell et al., 2004). The error propagation analysis of in situ well observations is even
more challenging due to temporal undersampling of the readings, and mischaracterization of aquifer properties. In this thesis, the error propagation scheme adopted is to express the uncertainty of a hypothetical, measured variable as its standard deviation so the final value of the variable and its associated error is defined as ( ). If variable follows a normal distribution, then one can be 68%
or ‘reasonably’ certain that the true value of lies within one (standard deviation) from the recorded value of . Similarly, one can be 95% or ‘strongly’ confident that the real value of lies in the region bounded by( ).
What follows is a short explanation of the principal uncertainty propagation techniques used in this work to compute uncertainties at each intermediate step all the way to the final results. First-order error propagation methods were used to compute uncertainties involving simple and straightforward calculations whereas Monte Carlo methods were resorted to in the case of more complex and nonlinear interaction of variables.
2.5.1. First-Order Error Propagation
Computation of uncertainty for simple arithmetic operations can be well approximated through first order error propagation formulae without recourse to the computationally intensive Monte Carlo simulations. Table 2.1 shows how to compute the standard deviation or uncertainty for simple operations. Here, represent the variables being operated on with standard deviations (errors) respectively. represent known scalar constants, and it is assumed that both are uncorrelated.
Table 2-1: First Order Error Propagation Approximation Formulae (Hiraishi et al., 2000)
Operation Uncertainty
√
√
√
2.5.2. Monte Carlo Method
Uncertainty analysis involving complex, nonlinear interaction between multiple variables cannot be well captured by first-order, linearized error propagation calculations. In this case, Monte Carlo analysis can be used to treat source variable uncertainties (Hiraishi et al., 2000). Simply put, Monte Carlo methods work in the following manner to provide bounds or distributions to an output function (Smith, 2006):
1) Choose a distribution that describes possible values of a parameter.
2) Generate synthetic data drawn randomly from this distribution.
3) Utilize the generated data as possible values of the input parameters to produce an output state space.
4) Study the distribution of the results in order to compute uncertainty. If the output is normally distributed or can be well approximated with a normal distribution, the mean value serves as the final output variable with the standard deviation as the uncertainty associated with the final result.
This study makes full use of Monte Carlo methods to generate synthetic GWS time series and explores the possible state space of both the intermediate and final results. For simple averaging and other operations, first-order error approximation techniques were instead implemented. In the case of Monte Carlo simulations, the information contained in the distribution of the output variables is then used to compute the resultant value and the uncertainty associated with it. More details about this procedure are contained in Chapter 4.
√(
) ( )
3. STUDY AREA
3.1. Background
The study area encompasses the Indian states predominantly located in the densely populated, heavily irrigated Indo-Gangetic plains. This region contains the states of Bihar, Haryana (including Delhi NCR), Punjab, Rajasthan, Uttar Pradesh, and West Bengal (see Fig 3-1). As of the Indian Census of 2011 (Chandramouli, 2011), the combined population of these states is around 530 million people and rising.
The state of Rajasthan indeed does not fall entirely into the Indo-Gangetic plains region except for its northern and eastern region but has been included in this study area as well. Lying between 68.2°-89.1° E Longitude, 24.3° - 32.2° N Latitude, the Indian portion of the Indo-Gangetic plains consists of the large floodplains of the Ravi, Beas, Sutlej, and Ganges Rivers, and are underlain by thick piles of Tertiary and Quaternary sediments (Jha & Sinha, 2009). This thick pile of alluvial deposits, which exceed 1000 meters in thickness at some locations, is extremely fertile and forms the largest consolidated area of irrigated food production in the world with a net cropped area of 114 million hectares (B. R. Sharma, Amarasinghe, &
Sikka, 2008). Indeed, more than 90% of the total water use in the area is for agriculture, with 8% followed by domestic use (B. Sharma et al., 2008).
Figure 3-1: Political map of India with the study area states labelled (“Political Map of India,” 2012)
Though the Indo-Gangetic plains region has a substantial surface-based irrigation infrastructure in place,
In fact, it is extremely difficult to find a farmer who either does not have his own pump or does not purchase water from the neighbouring pump owner (Jain & Singh, 2003). Groundwater pumping is one of the major driving forces behind the ‘green revolution’ of the 1970s which brought irrigation to areas not covered by surface irrigation, and helped catapult India from a nation reliant on food imports to a net food exporter (Shah, Roy, Qureshi, & Wang, 2003). The large-scale pumping of groundwater that fuelled the agricultural transformation has however led to increasing land and water degradation, water logging and salinization in highly irrigated areas, pollution of water resources, and declining water tables which all pose a grave threat to the water and food security of India.
3.2. Agriculture
Major crops grown in the study area are rice, wheat, cotton, millet, maize, and sugarcane which are grown in a dual cropping scheme: rice during the rabi period (October-March) and wheat during the kharif period (July-October) (Washington et al., 2012). Regional agricultural characteristics are detailed in Figure 3.2 with the western sub-regions 1, 2, and 3 (which roughly correspond to the states of North Rajasthan, Punjab, Haryana, Delhi NCR, Western U.P.) being a food surplus region featuring high productivity, high investment, and heavy use of groundwater for irrigation. In contrast, the eastern sub-regions 4, 5 (Eastern U.P., states of Bihar and West Bengal) form a largely food deficit region with low productivity, higher flood hazard risk, and poorer infrastructure (Aggarwal, Joshi, Ingram, & Gupta, 2004).
Figure 3-2: Agriculture zone map of Indo-Gangetic Plains (Aggarwal et al., 2004)
3.3. Climate
It is quite a challenge to generalize the climate of such a large study area however one can divide it into four seasons: Winter (December – April), Pre-Monsoon (April – June), Monsoon (June – September), and Post-Monsoon (October – December) (Sehgal, Singh, Chaudhary, Jain, & Pathak, 2013). The main source of rainfall in the region is the Southwest Monsoon which is known to account for 70-90% of the total annual rainfall over India (Rajeevan & McPhaden, 2004). Depending on the strength of the monsoon, typical annual rainfall over upper Indo-Gangetic Plains (Punjab, Haryana, North Rajasthan, and Western Uttar Pradesh) is 550 mm whereas lower Indo-Gangetic Plains (Eastern Uttar Pradesh, Bihar, West Bengal) is 1200 mm (Sehgal et al., 2013). However, this rainfall is unevenly distributed across the region with the dry Thar Desert in Rajasthan receiving less than 250 mm of rain in a year, and heavy rainfall being observed in parts of West Bengal to the tune of 2500 mm in a year (Central Ground Water Board, 2011).
3.4. Hydrogeology
The Indo-Gangetic plains are underlain with an extensive layer of Quaternary alluvial deposits, and bordered by the Himalayan mountain range to the north and the Deccan Shield to the south (Central Ground Water Board, 2011). The hydrogeological map of India is shown below in Figure 3-3. The region’s hydrogeology can be divided into three distinct, unconsolidated regions (Bhabar, Terai, Central Ganga Plains).
Bhabar Belt: Lying near the Himalayan foothills and confined to the northern parts of Punjab and Uttar Pradesh, the Bhabar Belt consists most of coarser alluvium (mainly river-deposited boulders and pebbles) forming the piedmont terrain (Central Ground Water Board, 2011). These coarse deposits are highly porous and allow for streams to flow easily underground thus allowing for extensive recharge of groundwater from the perennial rivers in the area. As a result, groundwater is unconfined and the water table is deep at 30 meters or beyond (Central Ground Water Board, 2011). The aquifers in this area are capable of yielding 100-300 m3/hour of water (Central Ground Water Board, 2011).
Terai Belt: Adjacent to the Bhabar Belt, the Terai Belt spans across Haryana, middle to upper Uttar Pradesh, Bihar, and upper West Bengal. The region consists of newer alluvium and is characterized by fine-grained sediments from the many rivers in the area (Central Ground Water Board, 2011). The geological structure of the Terai is highly porous and permeable thus the area has access to extensive groundwater resources. This results in an upper unconfined aquifer and a lower interconnected system of confined aquifers. Tubewells here report yields of 50-200 m3/hour of water (Central Ground Water Board, 2011).
Central Ganga Plain: Lying south of the Terai, the Central Ganga Plains contains perhaps the most productive aquifer system in India (Central Ground Water Board, 2011). Channel deposits and continuous flooding of the Ganges River creates a clastic, unconsolidated layer which hosts a widespread, multi- layered aquifer system underneath. The aquifer thickness varies from place to place, and range from a few meters to upwards of 300 meters. Though the aquifers are locally separated, they form a regional, hydraulically connected network (Karanth, 1987). Well yields here are typically in the 90 – 200 m3/hour range (Central Ground Water Board, 2011).
4. DATA AND METHODS
This thesis makes full use of remote sensing, land surface modelling, and ground-based observations to fulfil the research objectives. A short summary of the primary datasets and the CGWB-derived GWS solution is described in Table 4.1. The individual datasets and their respective pre-processing are detailed below.
4.1. GRACE-Derived Terrestrial Water Storage
As described in Section 2.3.1, GRACE Level-3 (1° x 1°, Monthly Resolution) Terrestrial Water Storage (TWS) Release-05 (RL-05) data product from the Jet Propulsion Laboratory (JPL) for January 2003 to December 2012 is used in this work. Level-3 GRACE products have all undergone major processing to correct for errors in certain spherical harmonics and post-glacial rebound. Spatial smoothing and destriping filters have further been applied to the data product (Landerer & Swenson, 2012; Swenson &
Wahr, 2006). The final output is expressed as TWS anomalies (in cm of equivalent water thickness) with respect to the average over January 2004 to December 2009. TWS anomalies during missing months have been gap filled using linear interpolation. The GRACE Level 3 solutions and the derivative data products listed below were all downloaded through GRACE TELLUS (“GRACE Land Mass Grids (Monthly),”
2013).
4.1.1. Land Grid Scaling
The previously mentioned filtering and de-striping operations are instrumental in reducing noise and removing certain correlated errors but also lead to loss in signal strength. To restore the geophysical signal, a 1° x 1° dimensionless, scaling coefficient grid is provided which needs to be multiplied with the corresponding land TWS grid (Landerer & Swenson, 2012). The scaled TWS anomaly is computed in the following manner:
( ) ( ) ( ) (4.1)
( ) represents each unscaled grid node, represents the longitude index, is the latitude index, is time (month) index, and ( ) is the scaling grid.
4.1.2. Error Handling
In addition to the land scaling grid, measurement and leakage error grids are also provided by GRACE TELLUS. These grids are also 1° by 1° in spatial resolution and are expressed in cm of equivalent water thickness. Measurement errors are due to inaccuracies inherent in the determination of the
