Repository containing scripts related to _Güntner et al., (2017, HESS)_ publication
======================================================================
**Article:**  
Andreas Güntner<sup>1,4</sup>, Marvin Reich<sup>1</sup>, Michal Mikolaj<sup>1</sup>, Benjamin Creutzfeldt<sup>2</sup>, Stephan Schroeder<sup>1</sup>, Hartmut Wziontek<sup>3</sup>: **_Landscape-scale water balance monitoring with an iGrav superconducting gravimeter in a field enclosure_**, published in: [HESS](http://www.hydrology-and-earth-system-sciences.net), **2017**

**Affiliation**  
<sup>1</sup>Helmholtz Centre Potsdam GFZ German Research Centre for Geosciences, Section 5.4 Hydrology, 14473 Potsdam, Germany  
<sup>2</sup>Senate Department for Urban Development and the Environment, 10707 Berlin, Germany  
<sup>3</sup>Federal Agency for Cartography and Geodesy (BKG), 04105 Leipzig, Germany  
<sup>4</sup>University of Potsdam, Institute of Earth and Environmental Science, 14476 Potsdam, Germany
> _Please cite this article when using here provided scripts_

## Dependencies
#### For Matlab scripts:
* Tested on Matlab R2015b + requires Statistics and Optimization Toolboxes
* Pull [igetsDataTools](https://github.com/emenems/igetsDataTools) to `\Scripts\MatlabLibrary\igetsDataTools`. _Tag:_ Guentner_et_al_2017
* Pull [hydroGravityLib](https://github.com/emenems/hydroGravityLib) to `\Scripts\MatlabLibrary\hydroGravityLib`. _Tag:_ Guentner_et_al_2017
* Pull [mGlobe](https://github.com/emenems/mGlobe) to `\Scripts\MatlabLibrary\mGlobe`

#### Hydrological model  
* Freely available [HYDRUS 1D](https://www.pc-progress.com/en/Default.aspx?hydrus-1d) was used for hydrological modelling. The executable in not provided here and must be [downloaded](https://www.pc-progress.com/en/Default.aspx?H1d-downloads)

#### Digital elevation model & GIS data
* Some scripts require local Digital elevation model or GIS data that are, however, not provided here due to provider's terms of use. Only sample DEM is provided here!

#### File formats
* _*.asc_ [arc ascii](https://en.wikipedia.org/wiki/Esri_grid) format   
* _*.tsf_ [TSoft](http://seismologie.oma.be/TSOFT/tsoft.html) file format
* _*.m_ scripts written in Matlab (R2015b). Some scripts may require Mapping, Statistics, Optimization Toolboxes.
* _*.r_ scripts written in R (3.3.1)

## Reproduction of results
####  Prepare gravity residuals
1. Register and Download [IGETS](http://isdc.gfz-potsdam.de/igets-data-base/) data. Preserve the data structure: `Data\Grav\IGETS\LevelX\YYYY` and file names!
	* Get iGrav-006 gravity time series: _Wettzell/we006/Level3/YYYY/IGETS-IGRAV-RESHOUR-we006-YYYYMMr1.ggp_
		* Download to `Data\Grav\IGETS\igrav006\Level3`
	* Get iGrav-006 auxiliary time series: _Wettzell/we006/Level3/YYYY/IGETS-IGRAV-AUX-we006-YYYYMMr1.ggp_
		* Download to `Data\Grav\IGETS\igrav006\Level3`
	* Get SG0XX time series: _Wettzell/we0XX-1/Level1/YYYY/IGETS-SG-MIN-we0XX-1-YYYYMM00.ggp_
		* Download to `Data\Grav\IGETS\sg0XX\Level1`
	* For the Noise level estimation, one second data are needed: _Wettzell/we0XX-1/Level1/YYYY/IGETS-XXX-SEC-we0XX-1-YYYYMM00.ggp_
		* Download to `Data\Grav\IGETS\xxxx\Level1\YYYY`
2. Compute iGrav006 gravity residuals:
	* Set and run `\Scripts\MatlabLibrary\igetsDataTools\igets_set_and_run.m` using `export_ggp = 1` option (only): go to _Export data to other formats_ and set
		* ```file_export = fullfile('..','..','..','\Data\Grav\','IGETS-IGRAV-RESHOUR-we006-ALL.tsf');```
		* ```path_export = fullfile('..','..','..','\Data\Grav\IGETS\igrav006\Level3\'); ```
		* ```prefix_export = 'IGETS-IGRAV-RESHOUR-we006-';```
		* ```suffix_export = 'r1.ggp';```
		* ```nanval_export = 99999.99;```
		* ```out_precision_export = '%10.2f';```
3. Compute SG030 gravity residuals
	* Set and run `\Scripts\MatlabLibrary\igetsDataTools\igets_set_and_run.m`
	* Converting 1 minute data to hourly residuals (`min_to_res = 1`): go to _One HOUR RESIDUAL data_ and set:
		* ```path_1hour = fullfile('..','..','..','Data','Grav','IGETS','sg030','Level3');```
		* ```prefix_1hour = 'IGETS-SG-RESHOUR-we030-1-';```
		* ```suffix_1hour = 'r1.ggp';```
		* ```prefix_1hour_log = 'IGETS-SG-RESHOUR-we030-1-';```
		* ```suffix_1hour_log = 'r1.log';```
		* ```calib_delay = -12.06;```
		* ```calib_factor = -739.519;```
		* ```drift_fit = [0.38796 -282868];```
		* ```input_tides = fullfile('..','..','..','\Data\Grav\','IGETS-IGRAV-RESHOUR-we006-ALL.tsf');```
		* ```input_tides_chan = 3```
		* ```tides_interp = 'spline';```
		* ```atmacs_admittance = -3;```
		* ```corr_file = fullfile('..','..','..','\Data\Grav\','SG030_correction_file_plotGrav.txt');```
		* ```fill_missing_1min = 60; ```
		* ```filter_file_1min = fullfile('data','g1m1h.nlf');```
		* atmacs_url_loc: [local link 1](http://atmacs.bkg.bund.de/data/results/lm/we_lm2_12km_19deg.grav) + [local link 2](http://atmacs.bkg.bund.de/data/results/iconeu/we_iconeu_70km.grav)
		* atmacs_url_glo: [global link 1](http://atmacs.bkg.bund.de/data/results/icon4lm/we_icon384_19deg.grav) + [global link 2](http://atmacs.bkg.bund.de/data/results/icon/we_icon384_20deg.grav)
	* Stacking monthly data to one file (`export_ggp = 1`): go to _Export data to other formats_ and set
		* ```file_export = fullfile('..','..','..','\Data\Grav\','IGETS-SG-RESHOUR-we030-1-ALL.tsf');```
		* ```path_export = fullfile('..','..','..','\Data\Grav\IGETS\sg030\Level3\'); ```
		* ```prefix_export = 'IGETS-SG-RESHOUR-we030-1-';```
		* ```suffix_export = 'r1.ggp';```
		* ```nanval_export = 99999.99;```
		* ```out_precision_export = '%10.2f';```
4. Combine iGrav006 and SG030 into one file and correct (iGrav) for PCB effect
	* Run `\Scripts\DataProcessing\igrav_correct_pcb.m`

####  Prepare other data
* Get hydro-meteorological time series: download from [IGETS](http://isdc.gfz-potsdam.de/igets-data-base/)
	* Get AUXiliary: _Wettzell/we006/Level3/YYYY/IGETS-IGRAV-AUX-we006-YYYYMMr1.aux_
		* Download to `Data\Grav\IGETS\igrav006\Level3\YYYY`
		* If in-situ meateorological measurements are available (not in IGETS), daily and hourly reference FAO ET can be computed using `\Scripts\MatlabLibrary\hydroGravityLib\et0calc.m`
	* Stacking monthly data to one file (`export_ggp = 1`): go to _Export data to other formats_ and set
		* ```file_export = fullfile('..','..','..','\Data\HydroMeteo\','IGETS-AUX-we006-ALL.tsf');```
		* ```path_export = fullfile('..','..','..','\Data\Grav\IGETS\igrav006\Level3\'); ```
		* ```prefix_export = 'IGETS-IGRAV-AUX-we006-';```
		* ```suffix_export = 'r1.aux';```
		* ```nanval_export = 99999.99;```
		* ```out_precision_export = '%10.4f';```
* To prepare data for optimization, run: `\Scripts\DataProcessing\Prepare_Optimization_Inputs.m`
* Get PCB temperature for May-August 2016
	* Download AUXiliary [IGETS](http://isdc.gfz-potsdam.de/igets-data-base/) data: _Wettzell/we006/Level1/YYYY/IGETS-IGRAV-AUX-we006-YYYYMM00.aux_
		* Download to `Data\Grav\IGETS\igrav006\Level1\2016`
	* Stacking monthly data to one file (`export_ggp = 1`): go to _Export data to other formats_ and set
		* ```file_export = fullfile('..','..','..','\Data\Grav\','iGrav006_PCB_temp_1min.tsf');```
		* ```path_export = fullfile('..','..','..','\Data\Grav\IGETS\igrav006\Level1\'); ```
		* ```prefix_export = 'IGETS-IGRAV-AUX-we006-';```
		* ```suffix_export = '00.aux';```
		* ```nanval_export = 99999.99;```
		* ```out_precision_export = '%10.4f';```
* Hydrological modeling (HYDRUS 1D)
    * all necessary input files to run the model are provided in `Data\HydrusModel\Input`
    * to run:
        * install [Hydrus 1D](https://www.pc-progress.com/en/Default.aspx?H1d-downloads)
        * create a new project
        * copy input files into this directory
        * run (with module or from within HYDRUS programm)
    * to re-create the synthetic input files run script `Scripts\HydrusModel\Hydrus_Inputdata.r`

#### Computations
* **Sensitivity to water storage variations**
	* Compute umbrella effect and surface area below building/pillar: Run `/Scripts/ForwardModel/Perform_Basic_Computations.m` set ```section.area``` and ```section.umbrella = 1``` (all other == 0)  
	* Integration radius effect (Figure8): Run `/Scripts/ForwardModel/Perform_Basic_Computations.m` set ```section.integ_initial``` and ```section.integ_sum = 1``` (all other == 0).
		* This script requires local Digital elevation model that is not provided here due to provider's conditions of use
		* Sample DEM only is provided here: `\Data\DEM\Sample_DEM.zip`
		* The results of the script used for publication are stored in `\Data\ForwardModel\All_IntegRadius_0007layers_radius1195m_umbella0200cm.mat`
		* Results will be used for Figure 8
	* Run virtual experiment using a hydrological model: see `Scripts/HydrusModel/` folder
	* Compute gravity effect per 1Vol%, as input for two artificial model runs and two gravimeter positions (iGrav006 and SG030): Run `/Script/ForwardModel/compute_gravity_response_per_1VolPerc.m`  
		* This script requires local Digital elevation model that is not provided here due to provider's conditions of use
		* The results of the script used for publication are stored in `\Data\ForwardModel\XXXX_gravity_response_1VolPerc_10m_deep_VertResolution_0p01m.txt`
		* Results will be used for Figure 9
	* To compute global hydrological (GHE) and non-tidal ocean loading effects (NTOL) pull and use [mGlobe](https://github.com/emenems/mGlobe)
		* All used GHE and NTOL time series are stored in this repository (`Data\Grav\GHE` & `Data\Grav\NTOL`)
* **Resolving water balance components**
	* Make sure all inputs have been created (see section above/Prepare other data)
	* For optimization of water budget parameters at annual scale run `Scripts\Optimization\Optimization_iGrav_TWS_annual.m`
	* For optimization of water budget parameters at daily scale run `Scripts\Optimization\daily_water_balance_solution.m`
	* To compute the RMSE of the daily scale approach, run `Scripts\Optimization\assess_optimization_performance_scale.m`

## Figures
* Figure 1: figure was prepared in QGIS. Input data are not provided here due to provider's terms of use
* Figure 2: open the `Scripts\Figures\Figure2.svg` file in inkscape (created using v0.92)
* Figure 3: run `Scripts\Figures\Figure3.m`
* Figure 4: run `Scripts\Figures\Figure4.m`
* Figures 5 & 6 were created by H. Wziontek, BKG, _hartmut.wziontek@bkg.bund.de_
* Figure 7: run `/Scripts/ForwardModel/Perform_Basic_Computations.m` while setting ```section.sensor_height = 1``` (all other == 0)
* Figure 8: run `Scripts\Figures\Figure8.m`
* Figure 9: run `Scripts\Figures\Figure9.m`
* Figure 10: run `Scripts\Optimization\Optimization_iGrav_TWS_annual.m`
* Figure 11: run `Scripts\Optimization\daily_water_balance_solution.m`
