HydroParam Preparation

1 Overview

This vignette documents the preparation and export of all static (time-invariant) input data for WaterGAP3. All outputs are written to the HydroParam/ directory as per-continent binary files (.fvec, .uvec, .umat, .field, .fmat), indexed by continent string idx_Conti.

2 CELL

Basic grid cell geometry at 1/12° resolution. Each cell receives a sequential Cell_ID following the order of entries in the file.

library(HydroGalleryCpp)
library(WaterGAP3Cpp)
library(WG3Tools)
library(tidyverse)

path_HydroParamExport <- "/projects/WaterGAP/WaterGAP3M/ProjectRun/Prepare_Iput/HydroParam/"

---

# CELL --------
## CELL_elevation_m ----------
map(str_Continent, \(idx_Conti) save_fvec(
  lst_CELL_elevation_m[[idx_Conti]],
  paste0(path_HydroParamExport, "CELL_elevation_m", "_", idx_Conti, ".fvec")
))

## CELL_latitude_Deg ---------
map(str_Continent, \(idx_Conti) save_fvec(
  lst_LonLat_WG3[[idx_Conti]]$y,
  paste0(path_HydroParamExport, "CELL_latitude_Deg", "_", idx_Conti, ".fvec")
))

## CELL_area_km2 --------
lst_CELL_area_km2 <- map(str_Continent, \(idx_Conti) grid_area_km2(lst_LonLat_WG3[[idx_Conti]]$x, lst_LonLat_WG3[[idx_Conti]]$y, num_Resolution_WG3))
map(str_Continent, \(idx_Conti) save_fvec(
  lst_CELL_area_km2[[idx_Conti]],
  paste0(path_HydroParamExport, "CELL_area_km2", "_", idx_Conti, ".fvec")
))

## CELL_aroundCell_ID ------------
# For each cell, stores the Cell_IDs of its 8 surrounding neighbours (Moore neighbourhood).
# Saved as an unsigned integer matrix (.umat).
lst_CELL_aroundCell_ID <- map(str_Continent, \(idx_Conti) grid_neighbour(lst_LonLat_WG3[[idx_Conti]]$x, lst_LonLat_WG3[[idx_Conti]]$y, num_Resolution_WG3))
map(str_Continent, \(idx_Conti) save_umat(
  lst_CELL_aroundCell_ID[[idx_Conti]],
  paste0(path_HydroParamExport, "CELL_aroundCell_ID", "_", idx_Conti, ".umat")
))

---

3 NET / RIVER

map(str_Continent, \(idx_Conti) save_uvec(
  lst_NET_downstreamCell_ID[[idx_Conti]],
  paste0(path_HydroParamExport, "NET_downstreamCell_ID", "_", idx_Conti, ".uvec")
))

Drainage network topology and upstream accumulation. The routing order is derived from the downstream cell map and stored as structured .field files rather than simple vectors, because the routing sequence is non-trivial (cells must be processed in upstream-to-downstream order).

# NET / RIVER -----------

## NET_routingStep_ID, NET_routingStepUpstream_ID -------------
# Derives the sequential routing order (step IDs) and the corresponding upstream cell
# assignments from the downstream cell map.
map(str_Continent, \(idx_Conti) generate_step_cell(lst_NET_downstreamCell_ID[[idx_Conti]],
                                                   paste0(path_HydroParamExport, "NET_routingStep_ID", "_", idx_Conti, ".field"),
                                                   paste0(path_HydroParamExport, "NET_routingStepUpstream_ID", "_", idx_Conti, ".field")))

## Upstream cells ----------
# Computes upstream cell lists and their total accumulated area (km²) per cell.
# Used downstream for land area correction and river length masking.
lst_upstreamCell_ID <- map(str_Continent, \(idx_Conti) get_inflow_cells(lst_NET_downstreamCell_ID[[idx_Conti]]))
lst_upstreamArea_km2 <- map(str_Continent, \(idx_Conti) {
  lst_Up_ID <- lst_upstreamCell_ID[[idx_Conti]]
  map_dbl(lst_Up_ID, function(idx_Id) {
    sum(lst_CELL_area_km2[[idx_Conti]][idx_Id], na.rm = TRUE)
  }) |> unlist()
})

---

4 Water Bodies

Five water body types are handled: Loclake (local lakes), Glolake (global lakes, directly connected to river routing), Locwetl (local wetlands), Glowetl (global wetlands), and Reservoi (reservoirs, human-managed, directly connected to river routing). A combined GlolakeReservoi object is also exported for the “Nature” run mode, where Glolakes and Reservoirs are treated identically.

4.1 Basic Area and Storage Capacity

Each water body type exports: cell_ID, area_km2, and capacity_m3.

walk(names(lst_Waterbodies_WG3), \(type) {

  ## cell_ID
  walk(str_Continent, \(idx_Conti) save_uvec(
    lst_Waterbodies_WG3[[type]][[idx_Conti]]$Cell_ID,
    paste0(path_HydroParamExport, type, "_cell_ID", "_", idx_Conti, ".uvec")
  ))

  ## area_km2
  walk(str_Continent, \(idx_Conti) save_fvec(
    lst_Waterbodies_WG3[[type]][[idx_Conti]]$Area_km2,
    paste0(path_HydroParamExport, type, "_area_km2", "_", idx_Conti, ".fvec")
  ))

  ## capacity_m3
  walk(str_Continent, \(idx_Conti) save_fvec(
    lst_Waterbodies_WG3[[type]][[idx_Conti]]$Capacity_m3,
    paste0(path_HydroParamExport, type, "_capacity_m3", "_", idx_Conti, ".fvec")
  ))
})

4.2 Reservoir Information

For each reservoir, two additional pieces of management information are required:

• dayControlStart_int: The day of the year (1–366) when reservoir management begins each year.
• isIrrigate_01: Indicates whether the reservoir is managed for irrigation (1) or not (0).

map(str_Continent, \(idx_Conti) save_uvec(
  as.integer(lst_Waterbodies_WG3$Reservoi[[idx_Conti]]$isIrrigate_01),
  paste0(path_HydroParamExport, "Reservoi_isIrrigate_01", "_", idx_Conti, ".uvec")
))

map(str_Continent, \(idx_Conti) save_uvec(
  lst_Waterbodies_WG3$Reservoi[[idx_Conti]]$dayControlStart_int,
  paste0(path_HydroParamExport, "Reservoi_dayControlStart_int", "_", idx_Conti, ".uvec")
))

The meanInflow_m3 can be obtained from either:

• the reservoir information provided in the input data, or
• the model-simulated natural runoff.

map(str_Continent, \(idx_Conti) save_fvec(
  lst_Waterbodies_WG3$Reservoi[[idx_Conti]]$meanInflow_m3,
  paste0(path_HydroParamExport, "Reservoi_meanInflow_m3", "_", idx_Conti, ".fvec")
))

Reservoir demand is allocated across up to 20 downstream grid cells. For each reservoir,

• downstream20Cell_ID gives the IDs of those cells, and
• downstream20AllocateFract_1 gives the fraction of demand assigned to each one.

map(str_Continent, \(idx_Conti) save_umat(
  lst_Reservoi_Downstream20_ID[[idx_Conti]],
  paste0(path_HydroParamExport, "Reservoi_downstream20Cell_ID", "_", idx_Conti, ".umat")
))

map(str_Continent, \(idx_Conti) save_fmat(
  lst_Reservoi_Downstream20_AllocateFract[[idx_Conti]],
  paste0(path_HydroParamExport, "Reservoi_downstream20AllocateFract_1", "_", idx_Conti, ".fmat")
))

4.3 Water Bodies Routing Step Fields

Each water body type also receives a routing step index field (.field), derived by intersecting the global routing order with the water body cell positions. These are used during time-stepping to process water bodies in correct topological order.

## Routing ------------
## NET_glolakeStepCell_ID
map(str_Continent, \(x) generate_step_extra_cell(paste0(path_HydroParamExport, "NET_routingStep_ID", "_", x, ".field"),
                                                 paste0(path_HydroParamExport, "Glolake_cell_ID", "_", x, ".uvec"),
                                                 paste0(path_HydroParamExport, "NET_glolakeStepCell_ID", "_", x, ".field")))
## NET_glowetlStepCell_ID
map(str_Continent, \(x) generate_step_extra_cell(paste0(path_HydroParamExport, "NET_routingStep_ID", "_", x, ".field"),
                                                 paste0(path_HydroParamExport, "Glowetl_cell_ID", "_", x, ".uvec"),
                                                 paste0(path_HydroParamExport, "NET_glowetlStepCell_ID", "_", x, ".field")))
## NET_reservoiStepCell_ID
map(str_Continent, \(x) generate_step_extra_cell(paste0(path_HydroParamExport, "NET_routingStep_ID", "_", x, ".field"),
                                                 paste0(path_HydroParamExport, "Reservoi_cell_ID", "_", x, ".uvec"),
                                                 paste0(path_HydroParamExport, "NET_reservoiStepCell_ID", "_", x, ".field")))
## NET_glolakeReservoiStepCell_ID
map(str_Continent, \(x) generate_step_extra_cell(paste0(path_HydroParamExport, "NET_routingStep_ID", "_", x, ".field"),
                                                 paste0(path_HydroParamExport, "GlolakeReservoi_cell_ID", "_", x, ".uvec"),
                                                 paste0(path_HydroParamExport, "NET_glolakeReservoiStepCell_ID", "_", x, ".field")))

---

5 Derived Cell Areas and River Length

After water body areas are known, land area and river length must be corrected. Land area subtracts the overlap of Glolakes/Reservoirs and Lakes from the full cell area. River length is set to zero for cells occupied by a Glolake or Reservoir, since those cells use a lake-type outflow scheme instead.

# Update Land area and river length ----------
## LAND_area_km2
lst_Land_Area <- map(str_Continent, \(idx_Conti) {
  
  vec_cell_area    <- lst_CELL_area_km2[[idx_Conti]]
  df_waterbody     <- lst_GlolakeReservoi_Area_WG3[[idx_Conti]]
  df_lake          <- lst_Loclake_WG3[[idx_Conti]]
  
  # Start with full cell area indexed by Cell_ID
  vec_land <- vec_cell_area
  
  # Subtract water body overlap area
  vec_land[df_waterbody$Cell_ID] <- vec_land[df_waterbody$Cell_ID] - df_waterbody$OverlapArea_km2
  
  # Subtract lake area
  vec_land[df_lake$Cell_ID] <- vec_land[df_lake$Cell_ID] - df_lake$Lake_area
  
  # Clamp to 0
  pmax(vec_land, 0)
})

map(str_Continent, \(idx_Conti) save_fvec(
  lst_Land_Area[[idx_Conti]],
  paste0(path_HydroParamExport, "LAND_area_km2", "_", idx_Conti, ".fvec")
))

## RIVER_length_km
## the cell with glolake or resevoi will set as zero
lst_RIVER_length_km_Korrekt <- map(str_Continent, \(idx_Conti) {
  river_Length <- lst_RIVER_length_km[[idx_Conti]]
  river_Length[lst_GlolakeReservoi_Area_WG3[[idx_Conti]]$Cell_ID] <- 0
  river_Length
})

map(str_Continent, \(idx_Conti) save_fvec(
  lst_RIVER_length_km_Korrekt[[idx_Conti]],
  paste0(path_HydroParamExport, "RIVER_length_km", "_", idx_Conti, ".fvec")
))

---

6 Physiographical Parameters

6.1 Climate Zone

Four parameters are spatially distributed by mapping climate zone class attributes onto the grid. Each parameter varies only across the discrete climate zone classes, not continuously.

# Climate zone ----------
## 4 Climate-zone-derived parameters for each continent -------------
str_Climate_Variable <- c(
  "param_EVATRANS_prt_alpha",
  "param_ATMOS_nettoRadiat_wg3a", "param_ATMOS_nettoRadiat_wg3b",
  "param_EVATRANS_wat_petmax"
)

map(str_Climate_Variable, \(idx_Vari) {
  lst_Data <- map_class_attribute(df_Climazone, lst_Climazone_WG3, idx_Vari)
  map(str_Continent, \(idx_Conti) save_fvec(
    lst_Data[[idx_Conti]],
    paste0(path_HydroParamExport, idx_Vari, "_", idx_Conti, ".fvec")
  ))
})

6.2 Land Cover

Seven parameters are distributed by land cover class. Root depth is also derived here and retained in memory for use in the soil capacity calculation below.

# Landcover ----------
## 7 landcover-class-derived parameters for each continent -------------
str_Landcover_Variable <- c(
  "LAND_builtRatio_1",
  "LAND_albedo_1", "LAND_snowAlbedo_1",
  "LAND_leafAreaIndexMin_1", "LAND_leafAreaIndexMax_1",
  "param_SNOW_fac_f"
)

map(str_Landcover_Variable, \(idx_Vari) {
  lst_Data <- map_class_attribute(df_Landcover_Hyde, lst_Landcover_WG3, idx_Vari)
  map(str_Continent, \(idx_Conti) save_fvec(
    lst_Data[[idx_Conti]],
    paste0(path_HydroParamExport, idx_Vari, "_", idx_Conti, ".fvec")
  ))
})

lst_LAND_growUpDay_d <- map_class_attribute(df_Landcover_Hyde, lst_Landcover_WG3, "LAND_growUpDay_d")
map(str_Continent, \(idx_Conti) save_uvec(
  lst_LAND_growUpDay_d[[idx_Conti]],
  paste0(path_HydroParamExport, "LAND_growUpDay_d", "_", idx_Conti, ".uvec")
))

## Root depth per landcover class, used to integrate soil water capacity by layer ----------
lst_RootDepth <- map_class_attribute(df_Landcover_Hyde, lst_Landcover_WG3, "LAND_rootDepth_m")

6.3 Soil

Three soil-derived fields are computed from the WISE30s database. All aggregations are weighted by soil component proportion (PROP) across layers.

# Soil ------------
## SOIL_capacity_mm: plant-available water capacity integrated over root depth -------
# weighted by soil component proportion (PROP) across WISE30s layers
df_Soil_Capacity <- df_Soil_WISE30s_DB |>
  select(SUID, SCID, PROP, Layer, TAWC) |>
  filter(TAWC > 0)

lst_WaterCapacity_WG3 <- soil_capacity_wise(df_Soil_Capacity, lst_Soil_WG3, lst_RootDepth) |>
  fill_na_wg3()

map(str_Continent, \(idx_Conti) save_fvec(
  lst_WaterCapacity_WG3[[idx_Conti]],
  paste0(path_HydroParamExport, "SOIL_capacity_mm", "_", idx_Conti, ".fvec")
))

## SOIL_potentialPercola_mm: maximum groundwater recharge rate -----------
# defined as the minimum bulk hydraulic conductivity across all soil layers
df_Soil_HydraulicConductivity <- df_Soil_WISE30s_DB |>
  filter(SDTO > 0, CLPC > 0, ORGC > 0, CFRAG > 0, BULK > 0) |>
  mutate(Data = hydrau_conduct_wise(SDTO, CLPC, ORGC, CFRAG, BULK)) |>
  select(SUID, SCID, PROP, Layer, Data)

lst_HydrauConduct_WG3 <- aggregate_soil_wise(df_Soil_HydraulicConductivity, lst_Soil_WG3, min) |>
  fill_na_wg3()

map(str_Continent, \(idx_Conti) save_fvec(
  lst_HydrauConduct_WG3[[idx_Conti]],
  paste0(path_HydroParamExport, "SOIL_potentialPercola_mm", "_", idx_Conti, ".fvec")
))

## USDA texture class per grid cell:  -----------
# weighted mean class index across soil components
# used to derive percolation parameters param_PERCOLA_wat_01 and param_PERCOLA_wat_k
df_Soil_ClassUSDA <- df_Soil_WISE30s_DB |>
  filter(SDTO > 0, CLPC > 0) |>
  mutate(Data = usda_texture_wise(SDTO, CLPC)) |>
  select(SUID, SCID, PROP, Layer, Data)

lst_ClassUSDA <- aggregate_soil_wise(df_Soil_ClassUSDA, lst_Soil_WG3, mean) |>
  fill_na_wg3()

## param_PERCOLA_wat_01: ------------
# binary flag, 1 if USDA class indicates low percolation (class index below threshold), 0 otherwise
lst_PERCOLA_wat_01 <- map(lst_ClassUSDA, \(x) as.integer(x < threshold_param_PERCOLA_wat_01_USDA))

map(str_Continent, \(idx_Conti) save_uvec(
  lst_PERCOLA_wat_01[[idx_Conti]],
  paste0(path_HydroParamExport, "param_PERCOLA_wat_01", "_", idx_Conti, ".uvec")
))

6.4 Slope — Percolation Scaling Factor

param_PERCOLA_wat_k is a cell-level percolation rate scaling factor assembled multiplicatively from four independent factor maps: USDA soil texture class, slope class, aquifer type, and permafrost extent.

# Slop -------
## param_PERCOLA_wat_k: -----------
#percolation rate scaling factor
# combined multiplicatively from soil texture, slope, aquifer, and permafrost factors
names(df_param_PERCOLA_wat_k_Factor_Soil)[1] <- "ID"
names(df_param_PERCOLA_wat_k_Factor_Slop)[1] <- "ID"

lst_Factor_Soil <- interpolate_class_attribute(df_param_PERCOLA_wat_k_Factor_Soil, lst_ClassUSDA,       "factor")
lst_Factor_Slop <- interpolate_class_attribute(df_param_PERCOLA_wat_k_Factor_Slop, lst_CELL_slop_perc, "factor")

lst_param_PERCOLA_wat_k <- pmap(
  list(lst_Factor_Soil, lst_Factor_Slop, lst_param_PERCOLA_wat_k_Aquifer, lst_param_PERCOLA_wat_k_Permafrost),
  \(a, b, d, e) a * b * d * e
)

map(str_Continent, \(idx_Conti) save_fvec(
  lst_param_PERCOLA_wat_k[[idx_Conti]],
  paste0(path_HydroParamExport, "param_PERCOLA_wat_k", "_", idx_Conti, ".fvec")
))

---

7 Global Constants (CONST)

Parameters in this section start as spatially uniform global constants and may be updated during calibration if identified as sensitive. They are broadcast to per-continent vectors of the appropriate length before export.

7.1 Clear-Sky Solar Radiation

Pre-computed for all 366 day-of-year values per cell, based on elevation and latitude. Stored as a float matrix (.fmat) with dimensions [n_cell × 366].

# CONST --------------

## ATMOS_solarRadiatClearSky_MJ ---------
map(str_Continent,
    \(idx_Conti) export_ATMOS_solarRadiatClearSky_MJ_366(paste0(path_HydroParamExport, "CELL_elevation_m", "_", idx_Conti, ".fvec"),
                                                 paste0(path_HydroParamExport, "CELL_latitude_Deg", "_", idx_Conti, ".fvec"),
                                                 paste0(path_HydroParamExport, "ATMOS_solarRadiatClearSky_MJ", "_", idx_Conti, ".fmat")))

7.2 Land CONST Parameters

The following parameters are set to a single global constant value and replicated across all cells per continent:

Parameter	Description
RIVER_velocity_km	Mean river flow velocity
param_INFILT_hbv_beta	HBV infiltration shape parameter
param_ATMOS_thr_Ts	Threshold temperature for snow/rain partitioning
param_SNOW_fac_Tmelt	Degree-day melt factor
param_EVATRANS_vic_gamma	VIC evapotranspiration parameter
param_EVATRANS_sup_k	Supplemental ET scaling factor
param_EVATRANS_sup_gamma	Supplemental ET shape parameter
param_PERCOLA_wat_thresh	Percolation threshold
param_BASEFLOW_sur_k	Surface baseflow recession coefficient

## LAND -----------
str_LAND_CONST_Variable <- c("RIVER_velocity_km",
                             "param_INFILT_hbv_beta",
                             "param_ATMOS_thr_Ts",
                             "param_SNOW_fac_Tmelt",
                             "param_EVATRANS_vic_gamma",
                             "param_EVATRANS_sup_k",
                             "param_EVATRANS_sup_gamma",
                             "param_PERCOLA_wat_thresh",
                             "param_BASEFLOW_sur_k")
map(str_LAND_CONST_Variable, 
    \(idx_Vari) {
      map(str_Continent, \(idx_Conti) save_fvec(
        rep(num_CONST_param_WG3[idx_Vari], n_Cell_Conti[idx_Conti]),
        paste0(path_HydroParamExport, idx_Vari, "_", idx_Conti, ".fvec")
      ))
      
    })

7.3 Water Body CONST Parameters

All three water body types share the same global albedo constant (Loclake_albedo_1). Outflow parameters differ by type: local water bodies use an area-capacity-power (acp) scheme, global water bodies use a linear storage (lin) scheme, wetland has also the snow process with kup scheme, and Reservoirs use the Hanasaki (han) scheme.

## Waterbody -----------
n_Loclake         <- map_int(str_Continent, \(idx_Conti) nrow(lst_Waterbodies_WG3$Loclake[[idx_Conti]]))
n_Glolake         <- map_int(str_Continent, \(idx_Conti) nrow(lst_Waterbodies_WG3$Glolake[[idx_Conti]]))
n_Locwetl         <- map_int(str_Continent, \(idx_Conti) nrow(lst_Waterbodies_WG3$Locwetl[[idx_Conti]]))
n_Glowetl         <- map_int(str_Continent, \(idx_Conti) nrow(lst_Waterbodies_WG3$Glowetl[[idx_Conti]]))
n_Reservoi        <- map_int(str_Continent, \(idx_Conti) nrow(lst_Waterbodies_WG3$Reservoi[[idx_Conti]]))
n_GlolakeReservoi <- map_int(str_Continent, \(idx_Conti) nrow(lst_Waterbodies_WG3$GlolakeReservoi[[idx_Conti]]))

save_const_param <- function(param_name, n_vec, out_name) {
  map(str_Continent, \(idx_Conti) save_fvec(
    rep(num_CONST_param_WG3[param_name], n_vec[idx_Conti]),
    paste0(path_HydroParamExport, out_name, "_", idx_Conti, ".fvec")
  ))
}

## Loclake
save_const_param("Water_albedo_1",                 n_Loclake, "Loclake_albedo_1")
save_const_param("param_Loclake_acp_storeFactor",   n_Loclake, "param_Loclake_acp_storeFactor")
save_const_param("param_Loclake_acp_gamma",         n_Loclake, "param_Loclake_acp_gamma")

## Glolake
save_const_param("Water_albedo_1",                 n_Glolake, "Glolake_albedo_1")
save_const_param("param_Glolake_lin_storeFactor",    n_Glolake, "param_Glolake_lin_storeFactor")

## Reservoi
save_const_param("Water_albedo_1",                 n_Reservoi, "Reservoi_albedo_1")
save_const_param("param_Reservoi_han_alpha",        n_Reservoi, "param_Reservoi_han_alpha")

## GlolakeReservoi
save_const_param("Water_albedo_1",                 n_GlolakeReservoi, "GlolakeReservoi_albedo_1")
save_const_param("param_Glolake_lin_storeFactor",    n_GlolakeReservoi, "param_GlolakeReservoi_lin_storeFactor")

## Locwetl
save_const_param("Water_albedo_1",                  n_Locwetl, "Locwetl_albedo_1")
save_const_param("param_Locwetl_acp_storeFactor",    n_Locwetl, "param_Locwetl_acp_storeFactor")
save_const_param("param_Locwetl_acp_gamma",          n_Locwetl, "param_Locwetl_acp_gamma")
save_const_param("param_Wetland_kup_Tfrozen",        n_Locwetl, "param_Locwetl_kup_Tfrozen")
save_const_param("param_Wetland_kup_f",              n_Locwetl, "param_Locwetl_kup_f")
save_const_param("param_Wetland_kup_Tmelt",          n_Locwetl, "param_Locwetl_kup_Tmelt")
save_const_param("param_Wetland_kup_MaxFrozenD",     n_Locwetl, "param_Locwetl_kup_MaxFrozenD")

## Glowetl
save_const_param("Water_albedo_1",                  n_Glowetl, "Glowetl_albedo_1")
save_const_param("param_Glowetl_lin_storeFactor",    n_Glowetl, "param_Glowetl_lin_storeFactor")
save_const_param("param_Wetland_kup_Tfrozen",        n_Glowetl, "param_Glowetl_kup_Tfrozen")
save_const_param("param_Wetland_kup_f",              n_Glowetl, "param_Glowetl_kup_f")
save_const_param("param_Wetland_kup_Tmelt",          n_Glowetl, "param_Glowetl_kup_Tmelt")
save_const_param("param_Wetland_kup_MaxFrozenD",     n_Glowetl, "param_Glowetl_kup_MaxFrozenD")