HAKESCH Calibration with new precipitation climatology from MeteoSwiss

Modifiziertes Fliesszeitverfahren

def modifiziertes_fliesszeitverfahren(
    x,              # Jährlichkeit: 2.3, 20, 100
    Vo20,           # Benetzungsvolumen 20-jährlich [mm]
    L,              # Fliesslänge [m]
    delta_H,        # Höhendifferenz entlang L [m]
    ys,             # Spitzenabflusskoeffizient [-]
    E,              # Einzugsgebietsfläche [km²]
    intensitaet_fn, # Funktion: i(x, Tc) in mm/h
    TB_start=10,    # Startwert für TB [min]
    tol=0.01,       # Konvergenztoleranz
    max_iter=100    # Max. Iterationen
):
    # 1. Benetzungsvolumen je nach x
    if x == 2.3:
        Vox = 0.5 * Vo20
    elif x == 100:
        Vox = 1.3 * Vo20
    elif x == 20:
        Vox = Vo20
    else:
        raise ValueError("Jährlichkeit x muss 2.3, 20 oder 100 sein.")

    # 2. Fliesszeit nach Kirpich
    J = delta_H / L
    TFl = 0.0195 * (L ** 0.77) * (J ** -0.385)

    # 3. Iteration zur Bestimmung von TB
    TB = TB_start
    for _ in range(max_iter):
        Tc = TB + TFl
        i = intensitaet_fn(x, Tc)  # [mm/h]
        TB_new = (Vox / i) * 60  # mm / (mm/h) -> h * 60 = min
        if abs(TB_new - TB) < tol:
            TB = TB_new
            break
        TB = TB_new
    else:
        raise RuntimeError("TB-Iteration konvergiert nicht.")

    Tc = TB + TFl
    i_final = intensitaet_fn(x, Tc)
    HQ = 0.278 * i_final * ys * E

    return {
        "HQ": HQ,
        "Tc": Tc,
        "TB": TB,
        "TFl": TFl,
        "i": i_final,
        "Vox": Vox
    }

---

Kölla-Verfahren

def koella_verfahren(
    x, LGe_km, Vo20, Gletscherflaeche_km2=0, rs=4,
    schneeschmelze=True, kGang=1.0, kFaktor=1.0,
    intensitaet_fn=None, TB_start=10, tol=0.01, max_iter=100
):
    FLeff = 0.12 * (LGe_km ** 0.17) * kFaktor  # in km²

    if x == 2.3:
        Vox = 0.5 * Vo20
        f = 0.1 * Vox
    elif x == 100:
        Vox = 1.3 * Vo20
        f = 0.1 * Vox
    elif x == 20:
        Vox = Vo20
        f = 0.1 * Vo20
    else:
        raise ValueError("Ungültige Jährlichkeit (x)")

    TFl_h = FLeff ** 0.2
    TFl = TFl_h * 60  # min

    TB = TB_start
    for _ in range(max_iter):
        Tc = TB + TFl
        i = intensitaet_fn(x, Tc) if intensitaet_fn else 100 * (Tc ** -0.4)
        TB_new = (Vox / i) * 60
        if abs(TB_new - TB) < tol:
            TB = TB_new
            break
        TB = TB_new
    else:
        raise RuntimeError("TB-Iteration konvergiert nicht.")

    Tc = TB + TFl
    i_final = intensitaet_fn(x, Tc) if intensitaet_fn else 100 * (Tc ** -0.4)
    if schneeschmelze:
        i_final += rs
    i_korrigiert = max(i_final - f, 0)

    QGle = 0.5 * Gletscherflaeche_km2

    HQ = FLeff * i_korrigiert * kGang + QGle

    return {
        "HQ": HQ,
        "Tc": Tc,
        "TB": TB,
        "TFl": TFl,
        "FLeff": FLeff,
        "i_final": i_final,
        "i_korrigiert": i_korrigiert,
        "QGle": QGle,
        "Vox": Vox,
        "Verlust": f
    }

---

#| label: Latin Hypercube Koella
# Load necessary libraries
library(reticulate)
library(lhs)        # For Latin Hypercube Sampling
library(ggplot2)    # For visualization

# Define target HQ value (e.g., from observations)
target_HQ <- 120.0

# Define the Latin-Hypercube Sampling (LHS) function to generate samples
lhs_samples <- randomLHS(10000, 2)  # Generate 10000 samples for kFaktor and kGang

# Rescale the LHS samples into the appropriate ranges for kFaktor and kGang
kFaktor_vals <- lhs_samples[, 1] * (3.0 - 0.1) + 0.1  # Range: 0.1 to 3.0
kGang_vals <- lhs_samples[, 2] * (2.0 - 0.5) + 0.5    # Range: 0.5 to 2.0

# Combine into a data frame for easier processing
samples_df <- data.frame(kFaktor = kFaktor_vals, kGang = kGang_vals)

# Run the function for each sample and compute the HQ
HQ_vals <- sapply(1:nrow(samples_df), function(i) {
  result <- py$koella_verfahren(
    x = 20,
    LGe_km = 50,
    Vo20 = 35,
    Gletscherflaeche_km2 = 1.2,
    rs = 4,
    schneeschmelze = TRUE,
    kGang = samples_df$kGang[i],
    kFaktor = samples_df$kFaktor[i]
  )
  return(result$HQ)  # Return HQ value from the function
})

# Add the computed HQ values to the data frame
samples_df$HQ <- HQ_vals

# Create a 2D scatter plot of kFaktor vs. kGang colored by HQ
ggplot(samples_df, aes(x = kFaktor, y = kGang, color = HQ)) +
  geom_point(alpha = 0.5) +
  scale_color_viridis_c() +  # Color scale for HQ
  labs(title = "HQ Values as Function of kFaktor and kGang",
       x = "kFaktor", y = "kGang", color = "HQ Value") +
  theme_minimal()

# Load necessary libraries
library(reticulate)
library(lhs)        # For Latin Hypercube Sampling
library(ggplot2)    # For visualization
library(viridis)    # For color scales

use_python(Sys.which("python"))  # Ensure Python environment is initialized
py_run_string("")  # Initialize Python access

# Define the Latin-Hypercube Sampling (LHS) function to generate samples
lhs_samples <- randomLHS(10000, 3)  # Generate 10000 samples for kFaktor, kGang, and Vo20

# Rescale the LHS samples into the appropriate ranges
kFaktor_vals <- lhs_samples[, 1] * (3.0 - 0.1) + 0.1   # kFaktor: Range 0.1 to 3.0
kGang_vals <- lhs_samples[, 2] * (2.0 - 0.5) + 0.5      # kGang: Range 0.5 to 2.0
Vo20_vals <- lhs_samples[, 3] * (100 - 10) + 10          # Vo20: Range 10 to 100

# Combine into a data frame for easier processing
samples_df <- data.frame(kFaktor = kFaktor_vals, kGang = kGang_vals, Vo20 = Vo20_vals)

# Run the function for each sample and compute the HQ
HQ_vals <- sapply(1:nrow(samples_df), function(i) {
  result <- py$koella_verfahren(
    x = 20,  # Fixed x value
    LGe_km = 50,
    Vo20 = samples_df$Vo20[i],
    Gletscherflaeche_km2 = 1.2,
    rs = 4,
    schneeschmelze = TRUE,
    kGang = samples_df$kGang[i],
    kFaktor = samples_df$kFaktor[i]
  )
  return(result$HQ)  # Return HQ value from the function
})

# Add the computed HQ values to the data frame
samples_df$HQ <- HQ_vals

# Now, we create 2D scatter plots for all pairs of parameters
# Function to create 2D scatter plot
plot_2d <- function(df, x_var, y_var) {
  ggplot(df, aes_string(x = x_var, y = y_var, color = "HQ")) +
    geom_point(alpha = 0.5) +
    scale_color_viridis_c() +  # Color scale for HQ
    labs(title = paste("HQ vs.", x_var, "and", y_var),
         x = x_var, y = y_var, color = "HQ Value") +
    theme_minimal()
}

# List of all parameter pairs to visualize
pairs_to_plot <- list(
  c("kFaktor", "kGang"),
  c("kFaktor", "Vo20"),
  c("kGang", "Vo20")
)

# Generate plots for each pair and display
for (pair in pairs_to_plot) {
  print(plot_2d(samples_df, pair[1], pair[2]))
}

---

📦 Use R to optimize kFaktor

library(reticulate)
library(ParBayesianOptimization)

use_python(Sys.which("python"))  # or use_virtualenv(), use_condaenv() if needed
py_run_string("")  # Initialize Python access

# Define objective function for calibration
obj_fun <- function(kFaktor) {
  result <- py$koella_verfahren(
    x = 20,
    LGe_km = 50,
    Vo20 = 35,
    Gletscherflaeche_km2 = 1.2,
    rs = 4,
    schneeschmelze = TRUE,
    kGang = 1.0,
    kFaktor = kFaktor
  )
  
  # Minimize the squared error to target HQ
  score <- (result$HQ - target_HQ)^2
  return(list(Score = score))
}

# Run Bayesian Optimization
opt_results <- bayesOpt(
  FUN = obj_fun,
  bounds = list(kFaktor = c(0.1, 3.0)),
  initPoints = 5,
  iters.n = 15,
  acq = "ucb"
)

opt_results$scoreSummary

opt_results$scoreSummary