PCR-GLOBWB 2 : a 5 arcmin global hydrological and water resources model

We present PCR-GLOBWB 2, a global hydrology and water resources model. Compared to previous versions of PCR-GLOBWB, this version fully integrates water use. Sector-specific water demand, groundwater and surface water withdrawal, water consumption, and return flows are dynamically calculated at every time step and interact directly with the simulated hydrology. PCR-GLOBWB 2 has been fully rewritten in Python and PCRaster Python and has a modular structure, allowing easier replacement, maintenance, and development of model components. PCR-GLOBWB 2 has been implemented at 5 arcmin resolution, but a version parameterized at 30 arcmin resolution is also available. Both versions are available as open-source codes on https://github. com/UU-Hydro/PCR-GLOBWB_model (Sutanudjaja et al., 2017a). PCR-GLOBWB 2 has its own routines for groundwater dynamics and surface water routing. These relatively simple routines can alternatively be replaced by dynamically coupling PCR-GLOBWB 2 to a global two-layer groundwater model and 1-D–2-D hydrodynamic models. Here, we describe the main components of the model, compare results of the 30 and 5 arcmin versions, and evaluate their model performance using Global Runoff Data Centre discharge data. Results show that model performance of the 5 arcmin version is notably better than that of the 30 arcmin version. Furthermore, we compare simulated time series of total water storage (TWS) of the 5 arcmin model with those observed with GRACE, showing similar negative trends in areas of prevalent groundwater depletion. Also, we find that simulated total water withdrawal matches reasonably well with reported water withdrawal from AQUASTAT, while water withdrawal by source and sector provide mixed results.


Introduction
The last decades saw the development of an increasing number of global hydrological models (GHMs), e.g. VIC (Liang et al., 1994;Nijssen et al., 2001), WMB (Fekete et al., 2002), WaterGAP (Döll et al., 2003), H08 (Hanasaki et al., Published by Copernicus Publications on behalf of the European Geosciences Union. 2430 E. H. Sutanudjaja et al.: PCR-GLOBWB 2 2008a, 2018, and Mac-PDM (Gosling and Arnell, 2011) (see Bierkens et al., 2014, and Kauffeldt et al., 2016, for a more extensive list, also including land surface models). GHMs have become essential tools to quantify and understand the global terrestrial water cycle, as they simulate the distributed hydrological response to weather and climate variations at higher resolution (typically 0.5 • × 0.5 • ) than used previously in general circulation models (GCMs), with more sophisticated run-off generation processes and river routing. As such, GHMs have been used for mediumrange to seasonal flood forecasting Alfieri et al., 2013;Candogan Yossef et al., 2013) as well as for a myriad of water-related global change assessments. Examples include the projection or estimation of future flood and drought events (Sperna-Weiland et al., 2012;Dankers et al., 2013;Prudhomme et al., 2013;Wanders et al., 2015;Wanders and Wada, 2016), current and future flood hazard and risk (Pappenberger et al., 2012;Hirabayashi et al., 2013;Winsemius et al., 2013Winsemius et al., , 2016, global groundwater depletion (Wada et al., 2010;Gleeson et al., 2012), the contribution of terrestrial water stores to global sea level change (Konikow, 2011;Pohkrel et al., 2013), current and future water scarcity under climate change and increasing population growth (Hanasaki et al., 2008b;Wada et al., 2011a, b;Schewe et al., 2014;Haddeland et al., 2014;, teleconnections between climate oscillations and water availability , the impact of land use change on global water resources (Rost et al., 2008;Sterling et al., 2015;Bosmans et al., 2017), and trends in surface water temperature and cooling water potential van Vliet et al., 2012). More recently, the output from global hydrological models has been extended to study socioeconomic impacts, such as virtual water trade (Konar et al., 2013;Dalin et al., 2017) and future agricultural production (Elliott et al., 2013). These applications show that GHMs have become invaluable tools in support of global change research and environmental assessments.
PCR-GLOBWB (PCRaster Global Water Balance)  is one of the recently developed GHMs. PCR-GLOBWB is a grid-based global hydrological model developed at the Department of Physical Geography, Faculty of Geosciences, Utrecht University, the Netherlands. The model, describing the terrestrial part of the hydrological cycle, was first introduced in a technical report by  and then formally published in a paper of , focusing on global water availability issues. PCR-GLOBWB was originally developed to solve the global daily surface water balance with a spatial resolution of 30 arcmin (about 50 km by 50 km at the Equator) and compare the resulting freshwater availability with monthly sectoral water demand in order to assess global-scale water scarcity Wada et al., 2011a, b). In this first version of PCR-GLOBWB (called PCR-GLOBWB 1 hereafter), similar to other global-scale hydrological models, water demand and water availability are treated independently, i.e. without direct feedback between human water use and other terrestrial water fluxes (e.g. Döll and Siebert, 2002;Wisser et al., 2010). Since it was first introduced, PCR-GLOBWB has been applied extensively in global water resource assessment studies. For instance, a recent search on Scopus (accessed on 13 April 2018) for the key-word "PCR-GLOBWB" yielded 113 publications with collectively over 2500 citations. Since the first version, several new model features have been introduced such as a comprehensive water demand and irrigation module (Wada et al., 2011b, a scheme for dynamic allocation of sectoral water demand to available surface water and groundwater resources, and the associated calculation of return flow . These features essentially introduced a two-way interaction among water demand, water withdrawal, water consumption, and availability, particularly over irrigated areas where water demand is large and return flow is significant. Nevertheless, all of these preceding studies using PCR-GLOBWB were performed at a relatively coarse resolution of 30 arcmin, limiting their subregional or local applications. Additionally, some added functionalities, such as the possibility to couple the land surface component of PCR-GLOBWB to a global MODFLOW-based groundwater model (Sutanudjaja et al., 2011de Graaf et al., , 2017 and an extension to simulate surface water temperature , were incorporated in different versions based on the original PCR-GLOBWB 1, leading to divergent model code development. The objective of this paper is to summarize and present the new version of the model, PCR-GLOBWB 2, which consolidates all components that have been developed since the original version of the model was first introduced . The new version of the model, PCR-GLOBWB 2, which is able to simulate the water balance at a finer spatial resolution of 5 arcmin, supersedes the original PCR-GLOBWB 1, which has a resolution of 30 arcmin only 1 . The finer resolution of PCR-GLOBWB 2 allows a much better representation of the effects of spatial heterogeneity in topography, soils, and vegetation on terrestrial hydrological dynamics (Wood et al., 2011;Bierkens et al., 2014). Likewise, it provides a better resolution for visualization that allows stakeholders and decision makers to assess model simulation output more easily and directly for the places they are specifically interested in (Sheffield et al., 2010;Beven and Cloke, 2012). To assess the possible improvements, this paper also presents the first evaluation results from the simulation of PCR-GLOBWB 2 at 5 arcmin resolution and compares them to a 30 arcmin version. As discharge data are commonly used in hydrological model performance evaluation, the simulated river discharge of PCR-GLOBWB 2 is compared to in situ discharge observations from the Global Runoff Data Centre (GRDC, 2014).
The paper is organized as follows. Section 2 provides a global description of PCR-GLOBWB 2, including its model structure and the new components and functionalities that have been added since PCR-GLOBWB 1. In Sect. 3 the global application of PCR-GLOBWB 2 is demonstrated and the results from a 58-year simulation  are evaluated against observations of discharge, total water storage, and reported withdrawal data. Section 4 summarizes and concludes this paper and discusses possible future developments. Section 5 provides information about availability of the model code and the underlying data.
2 PCR-GLOBWB 2 -model description 2.1 General overview PCR-GLOBWB 2 is a state-of-the-art grid-based global hydrology and water resources model. It is a component-based model implementation in Python using open-source PCRaster Python routines (Karssenberg et al., 2010, http://pcraster. geo.uu.nl/, last access: 15 September 2017. The code is distributed through GitHub. The computational grid covers all continents except Greenland and Antarctica. Currently two versions are available: one with a spatial resolution of 5 arcmin in latitude and longitude and one with a coarser resolution of 30 arcmin. Typical time steps for hydrology and water use are 1 day while sub-daily time stepping is used for hydrodynamic river routing. For all dynamic processes involved, PCR-GLOBWB 2 uses a time-explicit scheme. For each grid cell and each time step, PCR-GLOBWB 2 simulates moisture storage in two vertically stacked upper soil layers (S 1 + S 2 in Fig. 1), as well as the water exchange among the soil, the atmosphere, and the underlying groundwater reservoir (S 3 in Fig. 1). The exchange with the atmosphere is comprised of precipitation, evaporation from soils, open water, snow and soils, and plant transpiration, while the model also simulates snow accumulation and snowmelt. Sub-grid variability in land use, soils, and topography is included and influences the schemes for run-off-infiltration partitioning, interflow, groundwater recharge (from S 2 to S 3 ), and capillary rise (from S 3 to S 2 ). Run-off, generated by snowmelt, surface run-off, interflow, and baseflow, is routed across the river network to the ocean or endorheic lakes and wetlands. Routing can either be simple accumulation, simplified dynamic routing using a method of characteristics, or kinematic wave routing. In case the kinematic wave routing is used, it is also possible to use a (simplified) floodplain inundation scheme and to simulate the surface water temperature.
PCR-GLOBWB 2 includes a simple reservoir operation scheme that is applied to over roughly 6000 human-made reservoirs, which are progressively introduced according to Figure 1. Schematic overview of a PCR-GLOBWB 2 cell and its modelled states and fluxes. S 1 , S 2 (soil moisture storage), S 3 (groundwater storage), Q dr (surface run-off -from rainfall and snowmelt), Q sf (interflow or stormflow), Q bf (baseflow or groundwater discharge), and Inf (riverbed infiltration from to groundwater). The thin red lines indicate surface water withdrawal, the thin blue lines groundwater abstraction, the thin red dashed lines return flows from surface water use, and the thin dashed blue lines return flows from groundwater use surface. For each sector, withdrawal − return flow = consumption. Water consumption adds to total evaporation. In the figure, the five modules that make up PCR-GLOBWB 2 are portrayed on the model components.
their construction year, from the GRanD database . Human water use is fully integrated within the hydrological model, meaning that at each time step (1) water demands are estimated for irrigation, livestock, industry, and households, (2) these demands are translated into actual withdrawals from groundwater, surface water (rivers, lakes, and reservoirs), and desalinization, subject to availability of these resources and maximum groundwater pumping capacity in place, and (3) consumptive water use and return flows are calculated per sector.
As an option PCR-GLOBWB 2 can be partially or fully coupled to a two-layer global groundwater model based on MODFLOW (de Graaf et al., 2017). Recent work (Hoch et al., 2017a, b) also includes coupling PCR-GLOBWB 2 to either Delft3D Flexible Mesh (Kernkamp et al., 2011) or LISFLOOD-FP (Bates et al., 2010), which are model codes that can be used to solve the 1-D-2-D shallow water equations (or approximations thereof) for detailed inundation studies. PCR-GLOBWB 2 has a flexible modular structure (Fig. 1). The modular structure of PCR-GLOBWB 2, both in terms of model concepts and implementation (separate modules are called from a main program), makes it easy to modify or replace components according to specific objectives of the model application, to introduce new modules or components within the modelling system, and to couple it to existing codes.
There are currently five main hydrological modules in PCR-GLOBWB 2 as illustrated in Fig. 1 and briefly described in Sect. 2.3: meteorological forcing, land surface, groundwater, surface water routing, and irrigation and water use. For an extensive description of the underlying equations and methods used in each of these modules we refer to the following sources: the meteorological forcing module from van Beek (

Description of the modules
Hereafter, we briefly describe the main features of the five modules. Additionally, a (non-exhaustive) list of the model state and flux variables is provided in Table A1, whereas Table A2 lists the model inputs and parameters, including their sources.

Meteorological forcing module
Meteorological forcing of PCR-GLOBWB 2 uses time series of spatial fields of precipitation, temperature, and reference evaporation. Reference potential evaporation can be prescribed or calculated within the model and is used in the land surface module to calculate land-cover-specific potential evaporation based on crop factors of the various land cover types according to the FAO guidelines (Allen et al., 1998). There are two options for calculating reference potential evaporation: (1) using Hamon (1963) in case only daily mean temperature is available, or (2) using Penman-Monteith following the FAO guidelines (Allen et al., 1998) if net radiation, wind speed, and vapour pressure deficit are additionally available. See van Beek et al. (2008) for details. The resulting land-cover-specific potential evaporation is subsequently used to compute the actual evaporation for different land cover types in each cell. Apart from the calculation of evaporation, temperature is also used to partition precipitation into snow and rain and to drive snowmelt.

Land surface module
This core module of PCR-GLOBWB 2 covers the landatmosphere exchange, the vertical flow among soil compartments and the eventual groundwater recharge, snow and interception storage, and the run-off generation mechanisms. These processes are simulated over a number of land cover types and aggregated proportionally based on land cover fractions within a model cell. Users can specify their own land cover classification and introduce their own land cover parameterization. The number of land cover types is configurable. The standard parameterization of PCR-GLOBWB 2 carries four land cover types consisting of tall natural vegetation, short natural vegetation, non-paddy irrigated crops, and paddy irrigated crops (i.e. wet rice). There is also a parameterization set for six land cover types (Bosmans et al., 2017), albeit still at 30 arcmin resolution only, which includes distinct types for pasture and rain-fed crops. For the standard four land cover parameterization of PCR-GLOBWB, applied in this paper, the land cover types of pasture and rain-fed crops are integrated into the short natural vegetation type. For each land cover type, separate soil conditions can be specified. It should be noted that the soil and vegetation conditions are in any case fully spatially distributed. Thus, vegetation properties (e.g. crop factor; leaf area index, LAI) and soil properties (depth, saturated hydraulic conductivity, etc.) vary not only among land cover types but may also vary from cell to cell (e.g. per climate zone). In the standard parameterization, vegetation properties vary over the year using a monthly climatology of phenology and crop calendars (i.e. for the crop factor and LAI). The application of irrigation water for paddy and non-paddy irrigation is carried out by the irrigation and water use module. It is based on the FAO guidelines of Allen et al. (1998) and is dependent on the actual soil water storage (S 1 , S 2 ) or paddy-inundated water storages. All fluxes, from and to the land surface module in Fig. 1, are thus calculated separately per land cover type. The resulting vertical fluxes for each land cover type are interception evaporation, bare soil evaporation, snow sublimation, and vegetation-specific transpiration. In the soil column, vertical fluxes are driven by degrees of saturation of soil layers and interact with the underlying groundwater store, S 3 (see e.g. Sutanudjaja et al., 2011;Sutanudjaja, 2012, for detailed explanation). Surface run-off (Q dr , from precipitation and snowmelt) consists of infiltration excess run-off and saturation excess run-off following a sub-grid approach that mimics variable source areas, i.e. the improved Arno scheme (Todini, 1996;Hagemann and Gates, 2003). Interflow or stormflow (Q sf ), mostly occurring in regolith soils on hillslopes, is also handled with a sub-grid approach based on a run-off parameterization by Sloan and Moore (1984). All fluxes are computed per land cover type and balanced with the available storage to arrive at the net flux that is used to update the storages for the next time step. Also, to report the overall fluxes per cell, and to pass these to other modules, the land-cover-specific fluxes are subsequently averaged (weighted by land cover type fractions).
For the standard parameterization of the land surface module the following data sets are combined (see Table A2): the cell fractions of various non-irrigation land cover types are based on the map of Global Land Cover Characteristics Database (GLCC) version 2.0 (Loveland et al., 2000) with the land cover classification following Olson (1994a, b) and the parameter sets from  and Hagemann (2002). Irrigation land cover types (i.e. paddy and non-paddy), including their crop calendars and growing season lengths, are parameterized based on the data set of MIRCA2000  and the Global Crop Water Model of . We refer to  for detailed descriptions.

Groundwater module
The groundwater module calculates groundwater storage dynamics subject to recharge and capillary rise (calculated by the land surface module), groundwater discharge (Q bf , in case of a positive groundwater storage), and riverbed infiltration (Inf). Groundwater discharge (assumed to be the same as groundwater baseflow here) depends on a linear storageoutflow relationship (Q bf = S 3 /J ) in which the proportionality constant J is calculated following the drainage theory of Kraijenhoff van de Leur (1958) based on drainage network density and aquifer properties. Riverbed infiltration occurs only in the case that Q bf becomes 0 by groundwater withdrawal. Under persistent groundwater withdrawal (calculated with the irrigation and water use module) that is larger than the sum of recharge and riverbed infiltration, the groundwater storage S 3 is allowed to become negative. In this case, the part of the withdrawn groundwater in excess of the input (recharge and riverbed infiltration) is seen as non-renewable groundwater withdrawal leading to groundwater depletion (permanent loss of groundwater from storage). In case withdrawal becomes smaller than the input, the remaining input is used to first fill the negative storage to zero, before baseflow Q bf commences again. As an alternative, it is also possible to limit the maximum volume of non-renewable groundwater that can be extracted.
It is possible to use a full-fledged groundwater flow model based on MODFLOW (Harbaugh et al., 2000) coupled to PCR-GLOBWB 2 in order to calculate groundwater heads and flow paths. This can be performed as a one-way coupling in which PCR-GLOBWB 2 is first run with the standard groundwater module (reservoir S 3 with only vertical fluxes) to yield time series of net groundwater recharge (recharge -capillary rise) and surface water levels. These fluxes and inputs are subsequently used to force the groundwater flow model (see e.g. Sutanudjaja et al., 2011;de Graaf et al., 2017). Another possibility is to use a two-way coupling in which the groundwater module of PCR-GLOBWB 2 is replaced by the groundwater flow model. In this case, at each time step fluxes are exchanged between the groundwater model and the land surface module, and the groundwater model and the surface water routing module .

Surface water routing module
Following an eight-point steepest-gradient algorithm across the terrain surface (local drainage direction), all cells of the modelled domain are connected to a strictly convergent drainage network that together makes up the river basins and sub-basins of the model domain. The lowermost cell is either connected to the ocean or to an endorheic basin. Per cell, the sum of the three daily run-off fluxes ( Fig. 1) is aggregated and routed along the drainage network until passing the lowermost cell and being removed from the model. Routing can be carried out in three ways of increasing complexity: (1) simple accumulation of the fluxes over the drainage network, (2) a travel-time characteristic solution (Karssenberg et al., 2007), and (3) the kinematic wave solution.
The first method is typically aggregated over longer time steps (e.g. month or year) that are larger than the travel times of water along the longest river length. The second routing method includes an estimation of cell flow velocity based on average discharge from the last 5 years and Manning's equation, which assumes the energy slope to be equal to the bed slope. This estimated velocity is used to move the volume of water in the channel of a cell the corresponding distance within one daily time step along the drainage network. This method works reasonably well for relatively steep rivers in humid climates for which the friction slope is close to the bed slope and the rivers are equally filled with water throughout the year. The third method is the kinematic wave approximation of the Saint-Venant equations with flow described by Manning's equation. Also, here, it is assumed that friction slope and bed slope are equal, which makes it valid for rivers without backwater effects. The kinematic wave is solved using a time-explicit variable sub-time stepping scheme based on the minimum Courant number. Of these methods, the kinematic wave solution simulates the propagation of the flood wave more realistically while the others provide an expedient means to approximate discharge over longer periods.
Using the kinematic wave method, it is possible to model floodplain inundation, which occurs if the discharge exceeds the bankfull capacity of a channel. The excess discharge volume is spread over the entire cell from the lowest part of the cell (based on a higher-resolution sub-grid DEM) yielding a flooded area with an approximated flood depth. In case of flooding, the simulated river flow is impacted by adjusting the wetted area and wetted perimeter and calculating a weighted Manning coefficient from the individual Manning coefficients of the floodplains and the channel.
Lakes and reservoirs are part of the drainage network. Lakes and reservoirs can extend over multiple cells, in which case the storage is subdivided by area to ensure that lake and reservoir levels are the same across their extent. The active storage of lakes and the actual storage of reservoirs are dynamically updated; for the lake outflow a standard storageoutflow relationship based on a rectangular cross section over a broad-crested weir (Bos, 1989) is used, while reservoirs follow a release strategy. This strategy is, by default, aimed at passing the average discharge, while maintaining levels between a minimum and maximum storage , but more elaborate strategies that take account of downstream water demand are possible (e.g. . Lakes and reservoir areas change based on global volume-area relationships. All surface water areas, which can be classified into several water types including river channels, inundated floodplains, lakes and reservoirs, are subject to open water evaporation calculated from reference potential evaporation multiplied with factors depending on water types and depths. Moreover, surface waters are subject to surface water withdrawal calculated with the irrigation and water use module.
If the kinematic wave approach is used, it can also be augmented with an energy routing scheme to simulate surface water temperature . Finally, it should be noted that it is possible to run the routing routine from PCR-GLOBWB 2 as a stand-alone routine, which allows it to be fed with the specific discharge from other land surface models.
The routing methods that are available in PCR-GLOBWB 2 will yield significant errors for wide lowland rivers in which backwater effects are important. In this case, it is possible to replace the surface water module for part of the modelling domain with hydrodynamic models solving the shallow water equations .  developed a generic coupler for this purpose that enables coupling to multiple hydrodynamic modelling codes (https://doi.org/10.5281/zenodo.597107).
Although any data set can be used to define the drainage network and locate the lakes and reservoirs, the standard parameterization of PCR-GLOBWB 2 that runs globally uses the drainage network derived from the high-resolution 30 arcsec HydroSHEDS (Lehner et al., 2008) combined with 30 arcsec GTOPO30 (Gesch et al., 1999) and 1 km Hy-dro1k (Verdin and Greenlee, 1996;USGS EROS Data Center, 2006), lakes taken from the Global Lakes and Wetlands Database (GLWD) (Lehner and Döll, 2004) and reservoirs obtained from GRanD .

Irrigation and water use module
In PCR-GLOBWB 1 water demand was calculated separately from the hydrology and water availability calculated as a post-processing step by subtracting upstream demand (Wada et al., 2011a, b). In PCR-GLOBWB 2 water use (withdrawal and consumption) is fully integrated. Hereafter, the main features of the irrigation and water use module are described in the following order: water demand, water withdrawal, water consumption, and return flows.

Water demand
Irrigation water demand is calculated based on the crop composition (which changes per month and includes multicropping) and the irrigated area per cell. As stated above, these are obtained from MIRCA2000  and the Global Crop Water Model . In the standard PCR-GLOBWB 2 parameterization the irrigated areas change over time. In want of detailed data, fractions of paddy and non-paddy irrigation, as well as the crop composition per month, stay fixed (as obtained from MIRCA2000), while the total irrigated area per cell changes over time and is based on the FAOSTAT (http://www.fao.org/ faostat/en/\T1\textbackslash#home, last access: 15 September 2017) reported irrigated areas. Irrigation water demand is computed using the FAO guidelines (Doorenbos and Pruit, 1977;Allen et al., 1998): in case of non-paddy irrigation, water is applied whenever soil moisture falls below a pre-set value and then the soil column is replenished up to field capacity. In case of paddy irrigation, the water level is kept at a water depth of 5 cm above the surface until the late crop development stage (∼ 20 days) before the harvest. After that, no irrigation is applied anymore such that the water level is allowed to drop to zero under infiltration and evaporation . The net irrigation demand is augmented to account for limited irrigation efficiency and losses. In order to obtain irrigation water demand including losses, i.e. gross irrigation demand, net irrigation water demand is multiplied with (1 + f I ), with f I as a country-specific loss factor obtained from Rohwer et al. (2007).
Non-irrigation water demand covers three sectors: industry, households, and livestock. For each of these sectors, the gross demand and net demand are prescribed to the model. The calculation of net non-irrigation water demand, which varies with time, follows methods developed by . We refer to  for an extensive description. Trends in water demand are prescribed on an annual basis as a function of population, electricity demand, and gross domestic product (GDP) per capita. In addition, domestic water demand exhibits a seasonal variation on the basis of temperature. Domestic and industrial gross water demand is calculated from net water demand using a countryspecific recycling ratio (RC) (based on development stage or GDP per capita and additionally access to domestic water demand): gross = net / (1 − RC). This takes into account that much of the domestic and industrial water is not consumed but returned as surface water. For livestock, the return flow is assumed to be zero, meaning all water is consumed.

Water withdrawal
The water withdrawal estimation is based on the work by  and . In PCR-GLOBWB 2 water withdrawal is set equal to gross water demand (summed over all the sectors) unless sufficient water is not available. In that case, water withdrawal is scaled down to the available water and then allocated proportionally to gross water demand per sector. Thus, no allocation preference is available in the standard parameterization of PCR-GLOBWB 2.
Water can be abstracted from three sources: surface water, groundwater (fossil and non-fossil), and desalinated water. The latter is prescribed (Wada et al., 2011a), while the fractions of the other two sources are determined as a function of their relative abundance. Groundwater and surface water availability are determined based on 2-year running means of groundwater recharge and river discharge respectively, thus keeping track of the prevalence of local resources and their temporal change . These fractions determine, on a monthly basis, from which source water is abstracted. Surface water withdrawal is ceased if river discharge falls below 10 % of the long-term average yearly discharge under naturalized flow conditions (determined by running the model without withdrawal). If, for some reason, the surface water amount is insufficient, the model falls back on groundwater to meet the resulting gap. Groundwater is first abstracted from the renewable groundwater storage, and if this is not present, non-renewable groundwater is abstracted. The amount of groundwater that can be abstracted is, however, capped by the groundwater pumping capacity, which is based on data from the IGRAC GGIS database. The described dynamic allocation scheme is not always in line with local preferences or the infrastructure. However, there is a possibility to use fractions of groundwater and surface water withdrawal reported in the literature. For urban areas, we rely on the data set of McDonald et al. (2014) that states whether a surface water distribution infrastructure is available. If this is the case, industrial and domestic water withdrawals are mainly taken from surface water before abstracting groundwater. If surface water infrastructure is limited, groundwater source is prioritized (see e.g. Erkens and Sutanudjaja, 2015). For urban areas that are not in the McDonald (2014) data set, we give preference to the dynamic allocation scheme. For irrigation, we use the ratios supplied by  in regions where they are said to be reliable. In regions where they are not fully reliable, we take the average ratio provided by  and the one provided by the dynamic allocation scheme. For regions where the data of  are not reliable (i.e. extrapolated data), we give preference to the dynamic allocation scheme.
Moreover, we cannot assume that all the water demand is supplied from surface water and groundwater resources in the same cell. Ideally, data about local water redistribution networks and inter-basin transfers should be used to define surface water and groundwater service areas. Unfortunately, this information is not available at the global scale. Therefore, in our current parameterization of PCR-GLOBWB 2, we pool water availability of desalinated and surface water over zones of approximately 1 arcdeg by 1 arcdeg size that are truncated by country borders if applicable. For groundwater, 0.5 arcdeg zones are used. The downside of the current scheme is that a cell does not always have access to its nearest water resource if this lies outside its prescribed service area.

Water consumption and return flows
In case of irrigation, all the withdrawn water is applied to the soil (non-paddy) or the water level on the field (paddy). Part of that water is lost by transpiration and part by soil and open-water evaporation. Transpiration and evaporation together make up the irrigation water consumption. The remaining part of irrigated water is lost by percolation and contributes to groundwater recharge as return flow. Irrigation efficiency (not including conveyance losses) could also be calculated after the fact by the difference between withdrawal and transpiration. In the case of domestic and industrial water use, water consumption depends on the RC and equals withdrawal × (1 − RC), while withdrawal × RC constitutes return flow. All return flow is added to the surface water. For livestock, the consumption is set equal to the withdrawal and no return flow is assumed.

Model code
The original PCR-GLOBWB version 1  was written in the PCRaster scripting language. PCRaster (Wesseling et al., 1996) is a high-level programming language that started as a dynamic raster-based Geographical Information System (GIS) and is tailored to spa-tiotemporal modelling for environmental and earth science applications. The generic nature of PCRaster with its many pre-existing built-in hydrological functions and its syntax that reads like pseudo-code generally results in concise model codes, short development times, and limited programming errors. Karssenberg et al. (2010) developed a PCRaster Python package such that PCRaster functions, implemented in C++, can also be called via Python (http://www.python. org/). Using PCRaster Python makes it possible for students and beginner modellers to contribute to the model quickly, while it allows experts to be more productive and focus on the science rather than on the programming language syntax. Realizing the aforementioned advantages, PCR-GLOBWB, particularly starting from this version 2, has been rewritten in the Python scripting language.
To allow for exchanges of model components and therefore evaluate different model configurations, a componentbased development approach (e.g Argent, 2004;Castronova and Goodall, 2010) was followed while developing the PCR-GLOBWB 2 model code. Each of the PCR-GLOBWB scientific modules described in Sect. 2.3 is implemented in a separate Python class that needs to implement initialization and update methods. The latter designates changes of states and fluxes per time step. Each module is initialized and executed by iteratively calling the update method via a main model script.
To run the model, a so-called initialization file or configuration file is used (with extension .ini). In this file the following aspects are defined: the spatial and temporal domain, the time step, the settings of the different modules (e.g. with surface water routing, human water use, or not), and the locations and names of the parameter files and forcing files. PCR-GLOBWB 2 uses NetCDF files for most input and all output, thus making it easier to exchange data with other scientists and use existing tools to analyse their output.
PCR-GLOBWB 2 generally runs best under Linux. In order to run PCR-GLOBWB the following additional software needs to be installed: PCRaster version 4, Python version 2.7 with Python packages NumPy and netCDF4, and GDAL version 1.8 or higher.

Differences between PCR-GLOBWB 1 and 2
PCR-GLOBWB 2 has the following new capabilities compared to PCR-GLOBWB 1 (see Wada et al., 2011): the model was completely rewritten in PCRaster Python and now has a modular structure; the inputs and outputs are in the form of NetCDF files and output can be reported for daily monthly and yearly time steps; parameterizations are available at 30 and 5 arcmin resolutions; water use (demand, withdrawal, consumption, and return flow) is fully integrated; distinction is made between paddy and non-paddy irrigation and irrigation follows FAO guidelines; three different options for surface water routing are available and a surface water temperature module is fully integrated with the routing scheme; it is possible to run surface water routines separately with specific discharge from other sources (e.g. other land surface models); -PCR-GLOBWB 2 can be coupled to a two-layer transient groundwater model de Graaf et al., 2017) and to the hydrodynamic models Delft3D Flexible Mesh (Kernkamp et al., 2011) or LISFLOOD-FP (Bates et al., 2010;Hoch et al., 2017b).

Model demonstration and evaluation
To test and evaluate the performance of PCR-GLOBWB 2, we ran the model at both 30 and 5 arcmin resolution over the period 1958-2015. We compared the results of both simulations with discharge data from the Global Runoff Data Centre (GRDC, 2014), with total basin water storage estimates from GRACE (Gravity Recovery and Climate Experiment;  and with water withdrawal data from the FAO AQUASTAT database (FAO, 2016).

Parameterization
We used the standard parameterization (parameters, forcing, and their sources in Table A2) of PCR-GLOBWB 2 at 30 and 5 arcmin spatial resolutions to simulate global hydrology at daily resolution over 1958-2015. Outputs were reported as monthly averages. The parameterization was mostly unchanged from that given in van Beek and Bierkens (2009), but newer data sets, such as the GRAND  data set for reservoirs and MIRCA  for crop areas, were used if available. We stress that no calibration was performed. We ran the model with human water use options turned on and used the travel-time characteristic solution routing option.

Forcing
The forcing data set is based on time series of monthly precipitation, temperature, and reference evaporation from the CRU TS 3.2 data set of Harris et al. (2014) downscaled to daily values with ERA40 (1958-1978, Uppala et al., 2005 and ERA-Interim (1979, Dee et al., 2011. CRU is specified at 30 arcmin spatial resolution and directly usable. We used ERA40 and ERA-I results that had been resampled by ECMWF's resampling scheme from their original resolutions (∼ 1.2 and ∼ 0.7 • ) to 30 arcmin first. Here, resampling means a form of spatial downscaling whereby the values of the larger ERA40 and ERA-I grid cells are assigned to the cell centres and then spatially interpolated onto 30 arcmin grids. Precipitation was temporally downscaled by first applying a threshold of 0.1 mm day −1 to the ERA daily time series to estimate the number of rain days for ERA. The amount of rainfall below this threshold was proportionally allocated to the rain days. Next, the daily rainfall totals were scaled in order to reproduce the CRU monthly precipitation total using multiplicative scaling. Equally, monthly reference potential evaporation, computed with Penman-Monteith from the CRU data set, was scaled using multiplicative scaling and downscaled to daily data proportional to Hamon (1967) evaporation calculated from daily ERA temperatures. We elected not to calculate Penman-Monteith reference evaporation directly from the ERA40 and ERA-I data, in order to avoid the large calculation times needed to process the required meteorological values. For the air temperature, an additive scaling factor was used. To better simulate snow dynamics for the 5 arcmin model, the temperature values from CRU were further spatially downscaled to 5 arcmin using a temperature lapse rate derived from the higher-resolution CRU CL 2.0 climatology (New et al., 2002). For areas in which the number of stations underlying the CRU data set was found to be small, preference was given to directly using the meteorological data from ERA. The method used to create the forcing data set is described more extensively in van Beek (2008).

Spin-up
The large groundwater response times for certain regions (e.g. Niger and Amazon) requires substantial spin-up for the groundwater volumes to be in equilibrium with the current climate. To reach this equilibrium, the model was spun up using the average climatological forcing over the years 1958-2000 back to back for 150 years to reach a dynamic steady state. This spin-up was executed under naturalized conditions, which means no reservoirs and no human water use.

Computation time and parallelization
The models were run on Cartesius, the Dutch national supercomputer ( We used discharge stations from GRDC (2014) to compare simulated discharge from PCR-GLOBWB 2 with monthly reported discharge. From all the globally available stations in the database, we selected a subset of stations using the following criteria: (1) allowing a not-more-than 15 % difference in the catchment area between PCR-GLOBWB 2 and the area reported with the GRDC discharge station, (2) not more than 1 cell distance between the station location and the nearby location of a river in PCR-GLOBWB 2, and (3) at least 1 year of discharge data. This yielded 5363 stations for the 5 arcmin simulation, 3910 stations for the 30 arcmin simulation, and 3597 stations fulfilling the criteria for both resolutions. The minimum, median, and maximum catchment sizes for the GRDC stations at the 5 arcmin resolution are respectively 29, 2730, and 4.68 × 10 6 km 2 and 31, 6560, and 4.68 × 10 6 km 2 at the 30 arcmin resolution. As we jointly compared the performance of both simulations, we used the set of 3597 locations throughout. The average time series length of these stations is equal to 36 years.

Total water storage
We compared total water storage (TWS) as simulated by PCR-GLOBWB 2 with the TWS estimated from GRACE (Gravity Recovery and Climate Experiment) gravity anomalies. We used the JPL GRACE Mascon product RL05M Watkins et al., 2015;Wiese et al., 2016). Scanlon et al. (2016) suggest that recent developments in mascon (mass concentration) solutions for GRACE have significantly increased the spatial localization and amplitude of recovered terrestrial TWS signals. They also claim that one of the advantages of using the mascon solutions relative to traditional SH (spherical harmonic) solutions is that it makes it much easier for non-geodesists to apply GRACE data to hydrologic problems. Note that although the data of JPL RL05M are represented on a 30 arcmin lat-long grid, they represent the 3 × 3 arcdeg equal-area zones, which is the actual resolution of JPL RL05M. We compared trends on a pixel-bypixel basis. Given the coarse resolution of GRACE products of about 300 km by 300 km, we compared correlations only for major river basins with an area of 900 000 km 2 and up.

Water withdrawal
The water withdrawal for a large number of countries is taken from FAO's AQUASTAT database (FAO, 2016). These data are on average reported every 5 years. We compared simulated water withdrawal per sector and per water source (surface water and groundwater) with reported values per country and per reporting period, whenever available.

The global water balance simulated at 30 and 5 arcmin
We calculated the main global water balance components from the 30 and 5 arcmin simulations over the period 2000-2015. The results in Table 1 show that there are some differences between the two model runs, but values are on the same order of magnitude. The small difference in precipitation is due to the fact that the area of the land cells is slightly different at the two resolutions. Differences in evaporation and run-off show that the run-off and evaporation parameterization of PCR-GLOBWB 2 is not entirely scale consistent. Differences in evaporation may also be causing the differences in irrigation water demand, which in turn may explain the differences in water withdrawal. Recently, Samaniego et al. (2017) applied their multiscale parameter regionalization (MPR; creating spatially variable parameter fields) technique to PCR-GLOBWB 2 for the Rhine basin, showing that parameterizations that yield the same hydrological fluxes at different resolutions are possible. However, a global application of this method to all PCR-GLOBWB 2 parameters is not possible yet. Nonetheless, when comparing the results of both model runs with data reported in the literature, it shows that the global water balance components are similar to recent assessments (e.g. by Rodell et al., 2015), and groundwater withdrawal and total withdrawal estimates match those of previous studies (see Table 2). From Table 1, it can also be seen that there is a negative change in total terrestrial water storage in both model runs. Table 1 shows that this can only be partly explained by groundwater depletion, which is localized to certain regions (see also Sect. 3.4.2). Further analysis shows that this change can also be attributed to the trends in precipitation forcing used, particularly over the tropics.

Discharge
When evaluating the simulated discharge with discharge observations from GRDC, we used the monthly values and calculated three different measures. The first one is the correlation coefficient between monthly simulated and observed GRDC time series, which is a measure of reproducing correct timing of high and low discharge. A correlation coefficient of 1 indicates perfect timing. The second measure is the Kling-Gupta efficiency coefficient or KGE (Gupta et al., 2009), which equally measures bias, differences in amplitude, and differences in timing between monthly simulated and observed GRDC time series. The KGE varies between 1 and minus infinity, where 1 means a perfect fit in terms of bias, amplitude and timing. The last metric is the anomaly correlation, i.e. the correlation among monthly time series after the seasonal signal (climatology) has been removed. This statistic measures the ability of the model to correctly simulate timing of seasonal and inter-annual anomalies from the yearly climatology. This is to test if the model is able to capture the monthly-scale and inter-annual anomalies in discharge (i.e. on the monthly scale) when the dominant seasonal trend is removed from observations and simulations. An anomaly correlation of 1 indicates perfect characterization of inter-annual anomalies and values below 0 indicate a lack thereof. Figure 2 shows maps of the correlation coefficients for the GRDC stations considered and Fig. 3 shows histograms of correlation and KGE values. Both figures show that the evaluation results of the 5 arcmin simulation are generally better than those of the 30 arcmin simulation. For the 30 arcmin model, the number of catchments with KGE > 0, 0.3, and 0.6 are equal to 48, 26, and 7 % of the total catchments respectively. For the 5 arcmin model, these values are respectively equal to 63, 40, and 12 % of the total catchments. Note that for both runs the standard parameterization was used. Possible explanations for the better performance of the 5 arcmin run are a better delineation of the shape of the basins, particularly the smaller ones, a better characterization of basin relief and the drainage network, more accurate sub-grid parameterization of soil and land cover due to a smaller scale gap that needs to be overcome, better estimates of the basin storage, and better snow dynamics due to the downscaling of temperature to 5 arcmin resolution. The KGE values are less favourable than the correlation coefficients. This is mostly due to biases in run-off caused by incorrect meteorological forcing. It is difficult to exactly assess which of these factors are most important in determining the improvement. Inspecting the histograms of correlation and KGE (Fig. 3) shows that the improvement is mostly apparent for the smaller sized catchments, which supports the notion that a better delineation of the catchments' shape, topography, and drainage network could be the cause. However, disentangling these individual effects would require further study. To investigate the possible effects of better snow dynamics, we classified the GRDC stations into stations below 1000 m altitude (above mean sea level) and those above 1000 m. The GRDC stations above 1000 m are expected to experience precipitation falling as snow during periods of the year. The results in Fig. 4 clearly show that the improvement is larger for the higher GRDC stations. This supports the explanation that better snow dynamics due to temperature lapsing in combination with a better resolved digital elevation model is partly responsible for the superior results at 5 arcmin. We also investigated if improvements were notably different among climate zones, by separately calculating KGEs for GRDC stations in the Köppen-Geiger zones A (tropical), B (desert), C (temperate), and D (continental). The results (not shown) show that the improvement is equally visible for climate zones A, B, and C and less so for D (continental). Without further analysis this is difficult to explain. Note, however, that the continental climate zone is somewhat under-represented in the GRDC data set due to the low measurement densities over Russia, although it is well represented in the US. Thus, it may be that the global improvements shown in Fig. 3 are somewhat positively biased.
The maps of correlations (Fig. 2) show the best results in Europe and North America where the meteorological forcing is generally more accurate as a result of more data used in the reanalysis products and higher station availability in the CRU data set. Also, monsoon-dominated basins are well simulated  due to the strong seasonal nature of both forcing and related discharge. The improvement of the 5 arcmin simulation over the 30 arcmin simulation in Europe is mostly seen in the Alps and the Norwegian mountains. This reflects the fact that topography and thus snow dynamics is better represented at higher resolution as shown in Fig. 4. The least accurate results are obtained for some of the African rivers, in particular the Niger, where the groundwater recession coefficients are probably overestimated and inland delta evaporation is underestimated, for some rivers in the Rocky Mountains, which may be the result of errors in snow dynamics, and for continental eastern Europe, which is most likely explained by an overestimation of the groundwater recession constants.
The histograms of the anomaly correlation are shown in Fig. 5. The anomaly correlations are generally lower than the correlations, showing that seasonality explains part of the skill in many regions where seasonal variation is dominant when compared to intra-annual or inter-annual variability. Clearly, the 5 arcmin results are much better than those of the half-degree simulation, indicating a higher skill with regard to capturing inter-annual anomalies. Figure 6 shows a map of the difference between the anomaly correlation and the correlation for the 5 arcmin case. This map shows that there are some regions where the anomaly correlation is better than the correlation (blue colours), e.g. snow-dominated regions in Canada and the Niger basin. These are catchments where the model has difficulty reproducing the correct sea- Figure 6. Map showing for the 5 arcmin run the difference between the correlation and the anomaly correlation between simulated and observed discharge time series for 3597 GRDC discharge stations; negative values mean that the correlation is higher than the anomaly correlation.
sonality as a result of errors in snow dynamics (Canada) or groundwater dynamics (Niger). Also, in the case of the Niger River, not representing the inner delta flooding and resulting high evaporation may be the cause of poor seasonal timing of discharge. Generally, the PCR-GLOBWB 2 simulation is able to capture major groundwater-depleted regions as suggested by GRACE, such as those in the Central Valley aquifer, the High Plains aquifer, the North China Plain aquifer, and parts of the Middle East, Pakistan, and India. For these regions, the absolute rates of TWS change (i.e. TWS declines) of PCR-GLOBWB 2 are generally larger, while the spatial pattern in the GRACE map tends to be smoother. This is mainly due to the lower resolution and spatial averaging used in the GRACE product, as well as the fact that the current PCR-GLOBWB 2 simulation does not include lateral groundwater flow among cells. In the polar regions where GRACE estimates mass loss due to melting glaciers and ice sheets, PCR-GLOBWB 2 simulates accumulation as a result of a lack of glacier parameterization. Finally, there are some clear differences over the Amazon and some parts of Africa. A possible explanation are errors in meteorological forcing data, which are not very accurate in these parts, but also problems with the over-estimation of PCR-GLOBWB's groundwater response times in these regions, which therefore fail to be sufficiently sensitive to recent changes in terrestrial precipitation.

Total water storage
Further analyses were conducted at basin-scale resolution, for which both TWS time series of PCR-GLOBWB 2 and GRACE JPL RL05M were averaged over a river basin area map derived from the 5 arcmin PCR-GLOBWB drainage network. We identified all river basins with sizes larger than 900 000 km 2 , which is similar to the GRACE resolution. Smaller river basins were merged to the nearest river basins or grouped together. For the remaining map of large basins, the correlations between PCR-GLOBWB 2 and GRACE basin-average monthly and annual TWS time series were calculated. Monthly correlation provides information about PCR-GLOBWB's ability to correctly time TWS seasonal variability (with a value equal to 1 for perfect timing), while the correlation for annual time series measures inter-annual variability.
The results in Fig. 8 show that PCR-GLOBWB 2 is able to capture GRACE's TWS seasonality for most basins around the world, with the exception of some cold regions in high latitudes (e.g. the Yukon River basin, Iceland). This shortcoming is most likely due to the lack of a proper representation of glacier and ice processes in PCR-GLOBWB 2. As expected, the correlation values for inter-annual time series are generally lower than the ones for monthly time series. There are some areas with negative correlation values, such as the Amazon, Niger, and Nile river basins. Apart from the uncertainty in the GRACE signal, these deficiencies may be related to errors in model forcing and structural errors such as errors in the groundwater response time and the effects of wetlands that have not been represented sufficiently well.

Water withdrawal
We compared simulated water withdrawal data from PCR-GLOBWB 2 with reported withdrawal data per country from AQUASTAT (FAO, 2016). The results are shown subdivided per source (Fig. 9) and per sector (Fig. 10). Total water withdrawal and surface water withdrawal are simulated reasonably well (R 2 between 0.84 and 0.96 and regression slopes between 0.70 and 1.08). However, groundwater withdrawal is underestimated for the smaller water users. A likely explanation for this is occasional groundwater withdrawal by farmers during dry periods in areas that have not been mapped as irrigated crops in MIRCA, such as grasslands in Germany and the Netherlands, for example, while this groundwater withdrawal is reported in AQUASTAT.
When looking at water withdrawal per sector, results are mixed. The largest agricultural water users are well captured, but the smaller ones are clearly underestimated. This is related to the fact that in many regions of the smaller water use countries, water is used for irrigation only occasionally during dry summers, while these areas are not mapped as irrigated crops in MIRCA. Also, many of these countries use irrigation technology that is not part of MIRCA, e.g. subsurface drainage by artificially high surface water levels such as in a number of developed delta regions in the world. However, even though these smaller countries are not well represented, PCR-GLOBWB 2 is still able to capture the big water users, which have a significant impact on the water cycle and are most important for global-scale analyses.  Both industrial and domestic water withdrawals are underestimated. The underestimation of industrial water withdrawal is partly caused by the fact that we do not include water withdrawal for thermoelectric cooling of power plants. The underestimation of domestic water withdrawal comes from the fact that we assume that the priority of water allocation is proportional to demand. This means that in times of shortage, water withdrawal is reduced with an equal percentage for agriculture, industry, and domestic use. In many countries, however, there is a priority series, whereby domestic demand is first met, industrial demand next, and agricultural demand comes last. As a result, we underestimate domestic water withdrawal and it also partly causes the under-estimation of industrial water withdrawal. This is corroborated by plotting gross water demand (which would be withdrawal if no shortage would occur) against AQUASTAT data. These plots (not shown here) result in a regression slope of 0.68-0.75 for industrial demand and 0.78-0.92 for domestic demand. These results thus reveal that the water allocation scheme of PCR-GLOBWB 2 should be further improved.

Conclusions and future work
We presented the most recent version of the open-source global hydrology and water resource model PCR-GLOBWB. This version, PCR-GLOBWB 2, has a global coverage at  -1997 1998-2002 2003-2007 2008-2012 2013-2015 1993-1997 1998-2002 2003-2007 2008-2012 2013-2015  The uppermost plots (a, b) are for total water withdrawal, the middle ones (c, d) are groundwater withdrawal, and the lowermost charts (e, f) are surface water withdrawal. The regression coefficient is based on regression to non-log-transformed data with the intercept kept at zero.   -1997 1998-2002 2003-2007 2008-2012 2013-2015 R 2 1993-1997 1998-2002 2003-2007 2008-2012 2013-2015 R 2 1993-1997 1998-2002 2003-2007 2008- The uppermost plots (a, b) are for withdrawal for agricultural purposes, the middle ones (c, d) are industrial withdrawal, and the lowermost charts (e, f) are domestic. The regression coefficient is based on regression to non-log-transformed data with the intercept kept at zero. 5 arcmin resolution. Apart from the higher resolution, the new model has an integrated water use scheme, i.e. every day sector-specific water demand is calculated, resulting in groundwater and surface water withdrawal, water consumption, and return flows. Dams and reservoirs from the GRanD database  are added progressively according to their year of construction. PCR-GLOBWB 2 has been rewritten in Python and uses PCRaster Python functions (Karssenberg et al., 2007). It has a modular structure, which makes the replacement and maintenance of model parts easier. PCR-GLOBWB 2 can be dynamically coupled to a global two-layer groundwater model (de Graaf et al., 2017;Sutanudjaja et al., , 2011, and a one-way coupling to hydrodynamic models for large-scale inundation modelling  is also available.

1993
Comparing the 5 arcmin with 30 arcmin simulations using discharge data, we clearly find an improvement in the model performance of the higher-resolution model. We find a general increase in correlation, anomaly correlation, and KGE, indicating that the higher-resolution model is better able to capture the seasonality, inter-annual anomalies, and the general discharge characteristics. Also, PCR-GLOBWB 2 is able to reproduce trends and seasonality in total water storage as observed by GRACE for most river basins. It simulates the hotspots of groundwater decline that abound in GRACE as well. Simulated total water withdrawal matches reasonably well with reported water withdrawal from AQUASTAT, while water withdrawal by source and sector provide mixed results.
Future work will concentrate on further improving the water withdrawal and water allocation scheme, developing a full dynamic (two-way) coupling with hydrodynamic models, developing 5 and 1 km resolution (or higher) parameterizations of PCR-GLOBWB 2 using scale-consistent parameterizations (e.g. using MPR; Samaniego et al., 2017), incorporating a crop growth model, and solving the full surface energy balance. Other foreseeable developments are using the model in probabilistic settings and in data-assimilation frameworks.
Code and data availability. PCR-GLOBWB 2 is open source and distributed under the terms of the GNU General Public License version 3, or any later version, as published by the Free Software Foundation. The model code is provided through a GitHub repository: https://github.com/UU-Hydro/PCR-GLOBWB_model (Sutanudjaja et al., 2017a, https://doi.org/10.5281/zenodo.595656). This keeps users and developers immediately aware of any new revisions. Also, it allows developers to easily collaborate, as they can download a new version, make changes, and suggest and upload the newest revisions. The configuration INI files for the global 30 and 5 arcmin models and the associated model parameters and input files are provided at https://doi.org/10.5281/zenodo.1045338 (Sutanudjaja et al., 2017b). Development and maintenance of the official version (main branch) of PCR-GLOBWB 2 is conducted at the Department of Physical Geography, Utrecht University. Yet, contributions from external parties are welcome and encouraged. For news on the latest developments and papers published based on PCR-GLOBWB 2, we refer to http://www.globalhydrology.nl and for the underlying PCRaster Python code to http://pcraster.geo.uu. nl.
Appendix A Table A1. List (non-exhaustive) of state and flux variables defined in PCR-GLOBWB.