From 6bddd218081eb7a4dcef22656d0c3796f5d783db Mon Sep 17 00:00:00 2001 From: Dietmar Winkler Date: Tue, 31 Mar 2026 10:32:45 +0200 Subject: [PATCH 01/10] Combine the RunOff models Fixes #87 --- OpenHPL/Waterway/RunOff_zones.mo | 47 ++++++++++++++++++-------------- 1 file changed, 27 insertions(+), 20 deletions(-) diff --git a/OpenHPL/Waterway/RunOff_zones.mo b/OpenHPL/Waterway/RunOff_zones.mo index 3cbb226..13ea5ed 100644 --- a/OpenHPL/Waterway/RunOff_zones.mo +++ b/OpenHPL/Waterway/RunOff_zones.mo @@ -5,18 +5,18 @@ model RunOff_zones "Run off model. (with 10 height zones)" parameter Integer N = 10 "# of height zones" annotation ( Dialog(group = "Geometry")); // parameters of the hydrology model - parameter Modelica.Units.NonSI.Temperature_degC T_t=-3.91223331e-01 "Threshold temperature" annotation (Dialog(group="Physically-based parameters")); + parameter Modelica.Units.NonSI.Temperature_degC T_t=1 "Threshold temperature" annotation (Dialog(group="Physically-based parameters")); parameter SI.Area A[N] = ones(N) * 41.3e6 "Catchment area" annotation ( Dialog(group = "Geometry")); - parameter Real s_T = 3.99187122e-02 "Soil zone saturation threshold, m" annotation ( - Dialog(group = "Empirical parameters")), a_1 = 2.31870660e-06 "Discharge frequency for surface runoff, 1/sec" annotation ( - Dialog(group = "Empirical parameters")), a_2 = 4.62497942e-06 "Discharge frequency for fast runoff, 1/sec" annotation ( - Dialog(group = "Empirical parameters")), a_3 = 1.15749393e-06 "Discharge frequency for base runoff, 1/sec" annotation ( + parameter Real s_T = 20e-3 "Soil zone saturation threshold, m" annotation ( + Dialog(group = "Empirical parameters")), a_1 = 0.547 / 86400 "Discharge frequency for surface runoff, 1/sec" annotation ( + Dialog(group = "Empirical parameters")), a_2 = 0.489 / 86400 "Discharge frequency for fast runoff, 1/sec" annotation ( + Dialog(group = "Empirical parameters")), a_3 = 0.0462 / 86400 "Discharge frequency for base runoff, 1/sec" annotation ( Dialog(group = "Empirical parameters")), a_L[N] = {15.43, 3.97, 1.79, 0.81, 1.27, 1.44, 1.03, 2.32, 1.31, 0.57} .* 1e6 ./ A "Fractional area covered by lakes, -" annotation ( - Dialog(group = "Geometry")), g_T = 1.50394570e-01 "Ground saturation threshold, m" annotation ( - Dialog(group = "Physically-based parameters")), PERC = 5.79256691e-09 "preccolation from soil zone to base zone, m/sec" annotation ( - Dialog(group = "Empirical parameters")), beta = 1.00291469e+00 "Ground zone shape coefficient, -" annotation ( - Dialog(group = "Empirical parameters")), k_m = 3.47554015e-08 "Melting factor, m/deg/sec" annotation ( + Dialog(group = "Geometry")), g_T = 150e-3 "Ground saturation threshold, m" annotation ( + Dialog(group = "Physically-based parameters")), PERC = 0.6e-3 / 86400 "Percolation from soil zone to base zone, m/sec" annotation ( + Dialog(group = "Empirical parameters")), beta = 2 "Ground zone shape coefficient, -" annotation ( + Dialog(group = "Empirical parameters")), k_m = 4e-3 / 86400 "Melting factor, m/deg/sec" annotation ( Dialog(group = "Physically-based parameters")), PCORR = 1.05 "Precipitation correction - Rainfall, -" annotation ( Dialog(group = "Empirical parameters")), SCORR = 1.2 "Precipitation correction - Snowfall, -" annotation ( Dialog(group = "Empirical parameters")), CE = 0.04 "Model parameter for adjusted evapotranspiration, 1/deg" annotation ( @@ -42,6 +42,8 @@ model RunOff_zones "Run off model. (with 10 height zones)" Dialog(tab = "Evapotranspiration")), columns_month_temp[:] = 2:N + 1 "Columns with monthly average temperature variations for different height zones" annotation ( Dialog(tab = "Evapotranspiration")), columns_flow[:] = {2} "Column with real observed run off" annotation ( Dialog(tab = "Real run off")); + parameter Boolean useInput = false "If true, meteorological data is provided via input variables instead of data files" + annotation (Dialog(group = "Input"), choices(checkBox = true)); // variables SI.Height V_s_w[N] "Water content in soil zone", V_b_w[N] "Water content in base zone", V_g_w[N] "Water content in ground zone", V_s_d[N] "Dry snow"; //V_s_s[N] "Soggy snow"; @@ -52,11 +54,16 @@ model RunOff_zones "Run off model. (with 10 height zones)" SI.Velocity Vdot_p[N] "Precipitation"; Real a_e[N], a_sw[N], F_o, F_e, R2, err = 0.5e-3 "Small error, m"; // using data - Modelica.Blocks.Sources.CombiTimeTable temp_var(tableOnFile = true, columns = columns_temp, tableName = tableName_temp, fileName = fileName_temp); - Modelica.Blocks.Sources.CombiTimeTable prec_var(tableOnFile = true, fileName = fileName_prec, columns = columns_prec, tableName = tableName_prec); - Modelica.Blocks.Sources.CombiTimeTable evap_var(tableOnFile = true, fileName = fileName_evap, columns = columns_evap, tableName = tableName_evap); - Modelica.Blocks.Sources.CombiTimeTable month_temp(tableOnFile = true, fileName = fileName_month_temp, columns = columns_month_temp, tableName = tableName_month_temp); - Modelica.Blocks.Sources.CombiTimeTable flow_var(tableOnFile = true, fileName = fileName_flow, columns = columns_flow, tableName = tableName_flow); + Modelica.Blocks.Sources.CombiTimeTable temp_var(tableOnFile = true, columns = columns_temp, tableName = tableName_temp, fileName = fileName_temp) if not useInput; + Modelica.Blocks.Sources.CombiTimeTable prec_var(tableOnFile = true, fileName = fileName_prec, columns = columns_prec, tableName = tableName_prec) if not useInput; + Modelica.Blocks.Sources.CombiTimeTable evap_var(tableOnFile = true, fileName = fileName_evap, columns = columns_evap, tableName = tableName_evap) if not useInput; + Modelica.Blocks.Sources.CombiTimeTable month_temp(tableOnFile = true, fileName = fileName_month_temp, columns = columns_month_temp, tableName = tableName_month_temp) if not useInput; + Modelica.Blocks.Sources.CombiTimeTable flow_var(tableOnFile = true, fileName = fileName_flow, columns = columns_flow, tableName = tableName_flow) if not useInput; + input Real temp_input[N] if useInput "Zone temperature [degC]"; + input Real prec_input[N] if useInput "Zone precipitation [mm/day]"; + input Real evap_input if useInput "Potential evapotranspiration [mm/day]"; + input Real month_temp_input[N] if useInput "Monthly average zone temperature [degC]"; + input Real flow_input if useInput "Observed flow for R2 calculation [m3/s]"; // connector Modelica.Blocks.Interfaces.RealOutput Vdot_runoff annotation ( Placement(transformation(extent = {{90, -10}, {110, 10}}), iconTransformation(extent = {{80, -20}, {120, 20}}))); @@ -72,14 +79,14 @@ equation Vdot_tot = sum(A .* (Vdot_b2br + Vdot_s2sr + Vdot_s2fr)); for i in 1:N loop // Snow zone (Snow rourine) - T[i] = temp_var.y[i]; - Vdot_p[i] = prec_var.y[i] * 1e-3 / 86400; + T[i] = if useInput then temp_input[i] else temp_var.y[i]; + Vdot_p[i] = (if useInput then prec_input[i] else prec_var.y[i]) * 1e-3 / 86400; der(V_s_d[i]) = Vdot_p_s[i] - Vdot_d2w[i]; // + Vdot_w2d[i]; //der(V_s_s[i]) = Vdot_p_r[i] - Vdot_w2d[i] + Vdot_d2w[i] - Vdot_s2g[i]; Vdot_p_s[i] = if T[i] <= T_t then Vdot_p[i] * PCORR * SCORR * (1 - a_L[i]) else 0; Vdot_p_r[i] = if T[i] > T_t then Vdot_p[i] * PCORR * (1 - a_L[i]) else 0; - Vdot_d2w[i] = if T[i] > T_t and V_s_d[i] >= err * 1e-2 then k_m * (T[i] - T_t) * (1 - a_L[i]) else 0; + Vdot_d2w[i] = if T[i] > T_t and V_s_d[i] >= 0 then k_m * (T[i] - T_t) * (1 - a_L[i]) else 0; //Vdot_w2d[i] = if T[i]<=T_t and V_s_s[i]>=0 then k_m*(T_t-T[i]) else 0; //Vdot_s2g[i] = if T[i]>T_t and V_s_s[i]=err and V_s_d[i]>=err then (1+a_w)*Vdot_d2w[i]+Vdot_p_r[i] else 0; Vdot_s2g[i] = Vdot_p_r[i] + Vdot_d2w[i]; @@ -87,7 +94,7 @@ equation //der(V_g_w[i]) = Vdot_s2g[i] -Vdot_g2s[i] - (1-a)*Vdot_g_e[i]; der(V_g_w[i]) = Vdot_s2g[i] - Vdot_g2s[i] - a_e[i] .* Vdot_g_e[i]; Vdot_g2s[i] = if V_g_w[i] >= 0 and V_g_w[i] < g_T then (V_g_w[i] / g_T) ^ beta * Vdot_s2g[i] else Vdot_s2g[i]; - Vdot_epot[i] = evap_var.y[1] * 1e-3 / 86400 * (1 + CE * (T[i] - month_temp.y[i])); + Vdot_epot[i] = (if useInput then evap_input else evap_var.y[1]) * 1e-3 / 86400 * (1 + CE * (T[i] - (if useInput then month_temp_input[i] else month_temp.y[i]))); Vdot_g_e[i] = if V_g_w[i] < g_T then V_g_w[i] / g_T * Vdot_epot[i] else Vdot_epot[i]; a_e[i] = if V_s_d[i] < err then 1 else 0; // Soil zone (Upprec zone) @@ -103,8 +110,8 @@ equation Vdot_l_e[i] = a_L[i] * Vdot_epot[i]; end for; // Error - F_o = (flow_var.y[1] - 17.230144) ^ 2; - F_e = (flow_var.y[1] - Vdot_tot) ^ 2; + F_o = ((if useInput then flow_input else flow_var.y[1]) - 17.230144) ^ 2; + F_e = ((if useInput then flow_input else flow_var.y[1]) - Vdot_tot) ^ 2; R2 = 1 - F_e / F_o; Vdot_runoff = Vdot_tot; annotation (preferredView="info", From ae5915fa381a68e90644dde6e43f88f826573b0e Mon Sep 17 00:00:00 2001 From: Dietmar Winkler Date: Tue, 31 Mar 2026 12:02:35 +0200 Subject: [PATCH 02/10] Clean up of parameters and units --- OpenHPL/Waterway/RunOff_zones.mo | 195 +++++++++++++++++++------------ 1 file changed, 121 insertions(+), 74 deletions(-) diff --git a/OpenHPL/Waterway/RunOff_zones.mo b/OpenHPL/Waterway/RunOff_zones.mo index 13ea5ed..5af86e9 100644 --- a/OpenHPL/Waterway/RunOff_zones.mo +++ b/OpenHPL/Waterway/RunOff_zones.mo @@ -1,117 +1,164 @@ within OpenHPL.Waterway; model RunOff_zones "Run off model. (with 10 height zones)" extends OpenHPL.Icons.RunOff; - // height zone segmentation - parameter Integer N = 10 "# of height zones" annotation ( - Dialog(group = "Geometry")); - // parameters of the hydrology model - parameter Modelica.Units.NonSI.Temperature_degC T_t=1 "Threshold temperature" annotation (Dialog(group="Physically-based parameters")); - parameter SI.Area A[N] = ones(N) * 41.3e6 "Catchment area" annotation ( - Dialog(group = "Geometry")); - parameter Real s_T = 20e-3 "Soil zone saturation threshold, m" annotation ( - Dialog(group = "Empirical parameters")), a_1 = 0.547 / 86400 "Discharge frequency for surface runoff, 1/sec" annotation ( - Dialog(group = "Empirical parameters")), a_2 = 0.489 / 86400 "Discharge frequency for fast runoff, 1/sec" annotation ( - Dialog(group = "Empirical parameters")), a_3 = 0.0462 / 86400 "Discharge frequency for base runoff, 1/sec" annotation ( - Dialog(group = "Empirical parameters")), a_L[N] = {15.43, 3.97, 1.79, 0.81, 1.27, 1.44, 1.03, 2.32, 1.31, 0.57} .* 1e6 ./ A "Fractional area covered by lakes, -" annotation ( - Dialog(group = "Geometry")), g_T = 150e-3 "Ground saturation threshold, m" annotation ( - Dialog(group = "Physically-based parameters")), PERC = 0.6e-3 / 86400 "Percolation from soil zone to base zone, m/sec" annotation ( - Dialog(group = "Empirical parameters")), beta = 2 "Ground zone shape coefficient, -" annotation ( - Dialog(group = "Empirical parameters")), k_m = 4e-3 / 86400 "Melting factor, m/deg/sec" annotation ( - Dialog(group = "Physically-based parameters")), PCORR = 1.05 "Precipitation correction - Rainfall, -" annotation ( - Dialog(group = "Empirical parameters")), SCORR = 1.2 "Precipitation correction - Snowfall, -" annotation ( - Dialog(group = "Empirical parameters")), CE = 0.04 "Model parameter for adjusted evapotranspiration, 1/deg" annotation ( - Dialog(group = "Empirical parameters")); - //a_w = 0.03 "Saturation coeficieant, -" annotation (Dialog(group="Empirical parameters")), - //a = 0.001 "Snow surface fraction, -" annotation (Dialog(group="Empirical parameters")), - // parameters for outsource files with data - parameter String fileName_temp = Modelica.Utilities.Files.loadResource("modelica://OpenHPL/Resources/Tables/Temp_var.txt") "File with temperature variations in different height zones" annotation ( - Dialog(tab = "Temperature")), fileName_prec = Modelica.Utilities.Files.loadResource("modelica://OpenHPL/Resources/Tables/Prec_var.txt") "File with precipitation variations in different height zones" annotation ( - Dialog(tab = "Precipitation")), fileName_evap = Modelica.Utilities.Files.loadResource("modelica://OpenHPL/Resources/Tables/Evap_var.txt") "File with evapotranspiration variations during the year" annotation ( - Dialog(tab = "Evapotranspiration")), fileName_month_temp = Modelica.Utilities.Files.loadResource("modelica://OpenHPL/Resources/Tables/Month_av_temp.txt") "File with monthly average temperature for evapotranspiration calculation" annotation ( - Dialog(tab = "Evapotranspiration")), fileName_flow = Modelica.Utilities.Files.loadResource("modelica://OpenHPL/Resources/Tables/Flow_var_d.txt") "File with real observed run off" annotation ( - Dialog(tab = "Real run off")); - parameter String tableName_temp = "zones_temp" "Table with temperature variations in different height zones" annotation ( - Dialog(tab = "Temperature")), tableName_prec = "zones_prec" "Table with precipitation variations in different height zones" annotation ( - Dialog(tab = "Precipitation")), tableName_evap = "evap" "Table with evapotranspiration variations during the year" annotation ( - Dialog(tab = "Evapotranspiration")), tableName_month_temp = "month_temp" "Table with monthly average temperature" annotation ( - Dialog(tab = "Evapotranspiration")), tableName_flow = "flow" "Table with real observed run off" annotation ( - Dialog(tab = "Real run off")); - parameter Integer columns_temp[:] = 2:N + 1 "Columns with temperature variations for different height zones" annotation ( - Dialog(tab = "Temperature")), columns_prec[:] = 2:N + 1 "Columns with precipitation variations for different height zones" annotation ( - Dialog(tab = "Precipitation")), columns_evap[:] = {2} "Column with evapotranspiration variations during the year" annotation ( - Dialog(tab = "Evapotranspiration")), columns_month_temp[:] = 2:N + 1 "Columns with monthly average temperature variations for different height zones" annotation ( - Dialog(tab = "Evapotranspiration")), columns_flow[:] = {2} "Column with real observed run off" annotation ( - Dialog(tab = "Real run off")); + + parameter Integer N = 10 "Number of height zones" + annotation (Dialog(group = "Geometry")); + parameter SI.Area A[N] = ones(N) * 41.3e6 "Catchment area" + annotation (Dialog(group = "Geometry")); + parameter SI.PerUnit a_L[N] = {15.43, 3.97, 1.79, 0.81, 1.27, 1.44, 1.03, 2.32, 1.31, 0.57} .* 1e6 ./ A "Fractional area covered by lakes" + annotation (Dialog(group = "Geometry")); + + parameter Modelica.Units.NonSI.Temperature_degC T_t=1 "Threshold temperature" + annotation (Dialog(group="Physically-based parameters")); + parameter Real k_m(unit="m/deg/s") = 4e-3 / 86400 "Melting factor" + annotation (Dialog(group = "Physically-based parameters")); + parameter SI.Volume g_T = 150e-3 "Ground saturation threshold" + annotation (Dialog(group = "Physically-based parameters")); + + parameter SI.Volume s_T = 20e-3 "Soil zone saturation threshold" + annotation (Dialog(group = "Empirical parameters")); + parameter SI.Frequency a_1 = 0.547 / 86400 "Discharge frequency for surface runoff" + annotation (Dialog(group = "Empirical parameters")); + parameter SI.Frequency a_2 = 0.489 / 86400 "Discharge frequency for fast runoff" + annotation (Dialog(group = "Empirical parameters")); + parameter SI.Frequency a_3 = 0.0462 / 86400 "Discharge frequency for base runoff" + annotation (Dialog(group = "Empirical parameters")); + parameter SI.VolumeFlowRate PERC = 0.6e-3 / 86400 "Percolation from soil zone to base zone" + annotation (Dialog(group = "Empirical parameters")); + parameter Real beta = 2 "Ground zone shape coefficient" + annotation (Dialog(group = "Empirical parameters")); + parameter Real PCORR = 1.05 "Precipitation correction - Rainfall" + annotation (Dialog(group = "Empirical parameters")); + parameter Real SCORR = 1.2 "Precipitation correction - Snowfall" + annotation (Dialog(group = "Empirical parameters")); + parameter Real CE(unit="1/deg") = 0.04 "Model parameter for adjusted evapotranspiration" + annotation (Dialog(group = "Empirical parameters")); + parameter Boolean useInput = false "If true, meteorological data is provided via input variables instead of data files" - annotation (Dialog(group = "Input"), choices(checkBox = true)); - // variables - SI.Height V_s_w[N] "Water content in soil zone", V_b_w[N] "Water content in base zone", V_g_w[N] "Water content in ground zone", V_s_d[N] "Dry snow"; - //V_s_s[N] "Soggy snow"; + annotation (Dialog(tab = "Input data"), choices(checkBox = true)); + parameter String fileName_temp = Modelica.Utilities.Files.loadResource("modelica://OpenHPL/Resources/Tables/Temp_var.txt") "File with temperature variations in different height zones" + annotation (Dialog(tab = "Input data", group = "Temperature", enable=not useInput)); + parameter String tableName_temp = "zones_temp" "Table with temperature variations in different height zones" + annotation (Dialog(tab = "Input data", group = "Temperature", enable=not useInput)); + parameter Integer columns_temp[:] = 2:N + 1 "Columns with temperature variations for different height zones" + annotation (Dialog(tab = "Input data", group = "Temperature", enable=not useInput)); + + parameter String fileName_prec = Modelica.Utilities.Files.loadResource("modelica://OpenHPL/Resources/Tables/Prec_var.txt") "File with precipitation variations in different height zones" + annotation (Dialog(tab = "Input data", group = "Precipitation", enable=not useInput)); + parameter String tableName_prec = "zones_prec" "Table with precipitation variations in different height zones" + annotation (Dialog(tab = "Input data", group = "Precipitation", enable=not useInput)); + parameter Integer columns_prec[:] = 2:N + 1 "Columns with precipitation variations for different height zones" + annotation (Dialog(tab = "Input data", group = "Precipitation", enable=not useInput)); + + parameter String fileName_evap = Modelica.Utilities.Files.loadResource("modelica://OpenHPL/Resources/Tables/Evap_var.txt") "File with evapotranspiration variations during the year" + annotation (Dialog(tab = "Input data", group = "Evapotranspiration", enable=not useInput)); + parameter String tableName_evap = "evap" "Table with evapotranspiration variations during the year" + annotation (Dialog(tab = "Input data", group = "Evapotranspiration", enable=not useInput)); + parameter Integer columns_evap[:] = {2} "Column with evapotranspiration variations during the year" + annotation (Dialog(tab = "Input data", group = "Evapotranspiration", enable=not useInput)); + parameter String fileName_month_temp = Modelica.Utilities.Files.loadResource("modelica://OpenHPL/Resources/Tables/Month_av_temp.txt") "File with monthly average temperature for evapotranspiration calculation" + annotation (Dialog(tab = "Input data", group = "Evapotranspiration", enable=not useInput)); + parameter String tableName_month_temp = "month_temp" "Table with monthly average temperature" + annotation (Dialog(tab = "Input data", group = "Evapotranspiration", enable=not useInput)); + parameter Integer columns_month_temp[:] = 2:N + 1 "Columns with monthly average temperature variations for different height zones" + annotation (Dialog(tab = "Input data", group = "Evapotranspiration", enable=not useInput)); + + parameter String fileName_flow = Modelica.Utilities.Files.loadResource("modelica://OpenHPL/Resources/Tables/Flow_var_d.txt") "File with real observed run off" + annotation (Dialog(tab = "Input data", group = "Real run off", enable=not useInput)); + parameter String tableName_flow = "flow" "Table with real observed run off" + annotation (Dialog(tab = "Input data", group = "Real run off", enable=not useInput)); + parameter Integer columns_flow[:] = {2} "Column with real observed run off" + annotation (Dialog(tab = "Input data", group = "Real run off", enable=not useInput)); + + SI.Volume V_s_w[N] "Water content in soil zone"; + SI.Volume V_b_w[N] "Water content in base zone"; + SI.Volume V_g_w[N] "Water content in ground zone"; + SI.Volume V_s_d[N] "Dry snow"; SI.VolumeFlowRate Vdot_tot "Total runoff"; - SI.Velocity Vdot_s2b[N] "Runoff rate from soil zone to base zone", Vdot_pl[N] "Precipitation in lake", Vdot_b2br[N] "Runoff rate from base zone t obase runoff", Vdot_l_e[N] "Rate of evapotranspiration from lake", Vdot_g2s[N] "Runoff rate from ground zone to soil zone", Vdot_s2sr[N] "Runoff rate from soil zone to surface runoff", Vdot_s2fr[N] "Runoff rate from soil zone to fast runoff", Vdot_s2g[N] "Runoff rate from snow zone to ground zone", Vdot_g_e[N] "Evapotranspiration rate from ground zone", Vdot_p_r[N] "Precipitation in mainland in the form of snow", Vdot_d2w[N] "Melting rate from dry snow form to water snow form", Vdot_p_s[N] "Precipitation in mainland in the form of snow", Vdot_epot[N] "Evapotranspiration"; - //Vdot_w2d[N] "Freezing rate from water snow form to dry snow form"; + SI.VolumeFlowRate Vdot_s2b[N] "Runoff rate from soil zone to base zone"; + SI.VolumeFlowRate Vdot_pl[N] "Precipitation in lake"; + SI.VolumeFlowRate Vdot_b2br[N] "Runoff rate from base zone t obase runoff"; + SI.VolumeFlowRate Vdot_l_e[N] "Rate of evapotranspiration from lake"; + SI.VolumeFlowRate Vdot_g2s[N] "Runoff rate from ground zone to soil zone"; + SI.VolumeFlowRate Vdot_s2sr[N] "Runoff rate from soil zone to surface runoff"; + SI.VolumeFlowRate Vdot_s2fr[N] "Runoff rate from soil zone to fast runoff"; + SI.VolumeFlowRate Vdot_s2g[N] "Runoff rate from snow zone to ground zone"; + SI.VolumeFlowRate Vdot_g_e[N] "Evapotranspiration rate from ground zone"; + SI.VolumeFlowRate Vdot_p_r[N] "Precipitation in mainland in the form of snow"; + SI.VolumeFlowRate Vdot_d2w[N] "Melting rate from dry snow form to water snow form"; + SI.VolumeFlowRate Vdot_p_s[N] "Precipitation in mainland in the form of snow"; + SI.VolumeFlowRate Vdot_epot[N] "Evapotranspiration"; + Modelica.Units.NonSI.Temperature_degC T[N] "Ambient temperature"; - SI.Velocity Vdot_p[N] "Precipitation"; + SI.VolumeFlowRate Vdot_p[N] "Precipitation"; Real a_e[N], a_sw[N], F_o, F_e, R2, err = 0.5e-3 "Small error, m"; - // using data + Modelica.Blocks.Sources.CombiTimeTable temp_var(tableOnFile = true, columns = columns_temp, tableName = tableName_temp, fileName = fileName_temp) if not useInput; Modelica.Blocks.Sources.CombiTimeTable prec_var(tableOnFile = true, fileName = fileName_prec, columns = columns_prec, tableName = tableName_prec) if not useInput; Modelica.Blocks.Sources.CombiTimeTable evap_var(tableOnFile = true, fileName = fileName_evap, columns = columns_evap, tableName = tableName_evap) if not useInput; Modelica.Blocks.Sources.CombiTimeTable month_temp(tableOnFile = true, fileName = fileName_month_temp, columns = columns_month_temp, tableName = tableName_month_temp) if not useInput; Modelica.Blocks.Sources.CombiTimeTable flow_var(tableOnFile = true, fileName = fileName_flow, columns = columns_flow, tableName = tableName_flow) if not useInput; + input Real temp_input[N] if useInput "Zone temperature [degC]"; input Real prec_input[N] if useInput "Zone precipitation [mm/day]"; input Real evap_input if useInput "Potential evapotranspiration [mm/day]"; input Real month_temp_input[N] if useInput "Monthly average zone temperature [degC]"; input Real flow_input if useInput "Observed flow for R2 calculation [m3/s]"; - // connector - Modelica.Blocks.Interfaces.RealOutput Vdot_runoff annotation ( + + Modelica.Blocks.Interfaces.RealOutput Vdot_runoff "Output connector" + annotation ( Placement(transformation(extent = {{90, -10}, {110, 10}}), iconTransformation(extent = {{80, -20}, {120, 20}}))); initial equation - //V_s_s = zeros(N); V_s_d = zeros(N); V_g_w = zeros(N); V_s_w = zeros(N); V_b_w = zeros(N); equation - // Total runoff - //der(Vdot_tot) = if V_s_w > s_T then (a_1*(V_s_w-s_T)+a_2*V_s_w)*(1-sum(a_L)/N)*sum(A) + sum(A)*a_3*V_b_w else a_2*V_s_w*(1-sum(a_L)/N)*sum(A) + sum(A)*a_3*V_b_w; - Vdot_tot = sum(A .* (Vdot_b2br + Vdot_s2sr + Vdot_s2fr)); + Vdot_tot = sum(A .* (Vdot_b2br + Vdot_s2sr + Vdot_s2fr)) "Total runoff"; for i in 1:N loop - // Snow zone (Snow rourine) - T[i] = if useInput then temp_input[i] else temp_var.y[i]; - Vdot_p[i] = (if useInput then prec_input[i] else prec_var.y[i]) * 1e-3 / 86400; + if useInput then + T[i] = temp_input[i]; + Vdot_p[i] = prec_input[i] * 1e-3 / 86400; + else + T[i] = temp_var.y[i]; + Vdot_p[i] = prec_var.y[i] * 1e-3 / 86400; + end if; der(V_s_d[i]) = Vdot_p_s[i] - Vdot_d2w[i]; - // + Vdot_w2d[i]; - //der(V_s_s[i]) = Vdot_p_r[i] - Vdot_w2d[i] + Vdot_d2w[i] - Vdot_s2g[i]; Vdot_p_s[i] = if T[i] <= T_t then Vdot_p[i] * PCORR * SCORR * (1 - a_L[i]) else 0; Vdot_p_r[i] = if T[i] > T_t then Vdot_p[i] * PCORR * (1 - a_L[i]) else 0; Vdot_d2w[i] = if T[i] > T_t and V_s_d[i] >= 0 then k_m * (T[i] - T_t) * (1 - a_L[i]) else 0; - //Vdot_w2d[i] = if T[i]<=T_t and V_s_s[i]>=0 then k_m*(T_t-T[i]) else 0; - //Vdot_s2g[i] = if T[i]>T_t and V_s_s[i]=err and V_s_d[i]>=err then (1+a_w)*Vdot_d2w[i]+Vdot_p_r[i] else 0; Vdot_s2g[i] = Vdot_p_r[i] + Vdot_d2w[i]; - // Ground zone (Soil moisure) - //der(V_g_w[i]) = Vdot_s2g[i] -Vdot_g2s[i] - (1-a)*Vdot_g_e[i]; - der(V_g_w[i]) = Vdot_s2g[i] - Vdot_g2s[i] - a_e[i] .* Vdot_g_e[i]; + + der(V_g_w[i]) = Vdot_s2g[i] - Vdot_g2s[i] - a_e[i] .* Vdot_g_e[i] "Ground zone (Soil moisure)"; Vdot_g2s[i] = if V_g_w[i] >= 0 and V_g_w[i] < g_T then (V_g_w[i] / g_T) ^ beta * Vdot_s2g[i] else Vdot_s2g[i]; - Vdot_epot[i] = (if useInput then evap_input else evap_var.y[1]) * 1e-3 / 86400 * (1 + CE * (T[i] - (if useInput then month_temp_input[i] else month_temp.y[i]))); + if useInput then + Vdot_epot[i] = evap_input * 1e-3 / 86400 * (1 + CE * (T[i] - month_temp_input[i])); + else + Vdot_epot[i] = evap_var.y[1] * 1e-3 / 86400 * (1 + CE * (T[i] - month_temp.y[i])); + end if; Vdot_g_e[i] = if V_g_w[i] < g_T then V_g_w[i] / g_T * Vdot_epot[i] else Vdot_epot[i]; a_e[i] = if V_s_d[i] < err then 1 else 0; - // Soil zone (Upprec zone) - der(V_s_w[i]) = Vdot_g2s[i] - a_sw[i] * Vdot_s2b[i] - Vdot_s2sr[i] - Vdot_s2fr[i]; + + der(V_s_w[i]) = Vdot_g2s[i] - a_sw[i] * Vdot_s2b[i] - Vdot_s2sr[i] - Vdot_s2fr[i] "Soil zone (Upprec zone)"; Vdot_s2b[i] = (1 - a_L[i]) * PERC; Vdot_s2sr[i] = if V_s_w[i] > s_T then a_1 * (V_s_w[i] - s_T) else 0; Vdot_s2fr[i] = a_2 * V_s_w[i]; a_sw[i] = if Vdot_g2s[i] < Vdot_s2b[i] then 0 else 1; - // Basement zone (Lower zone) - der(V_b_w[i]) = Vdot_s2b[i] + Vdot_pl[i] - Vdot_b2br[i] - Vdot_l_e[i]; + + der(V_b_w[i]) = Vdot_s2b[i] + Vdot_pl[i] - Vdot_b2br[i] - Vdot_l_e[i] "Basement zone (Lower zone)"; Vdot_pl[i] = a_L[i] * Vdot_p[i]; Vdot_b2br[i] = a_3 * V_b_w[i]; Vdot_l_e[i] = a_L[i] * Vdot_epot[i]; end for; - // Error - F_o = ((if useInput then flow_input else flow_var.y[1]) - 17.230144) ^ 2; - F_e = ((if useInput then flow_input else flow_var.y[1]) - Vdot_tot) ^ 2; + + if useInput then + F_o = (flow_input - 17.230144) ^ 2; + F_e = (flow_input - Vdot_tot) ^ 2; + else + F_o = (flow_var.y[1] - 17.230144) ^ 2; + F_e = (flow_var.y[1] - Vdot_tot) ^ 2; + end if; R2 = 1 - F_e / F_o; Vdot_runoff = Vdot_tot; annotation (preferredView="info", @@ -324,4 +371,4 @@ recession constant for the base runoff in lower zone. Finally, the user can also the data for the CombiTimeTable models are stored.

")); -end RunOff_zones; +end RunOff_zones; \ No newline at end of file From b7501b8a0bb73eaf7a0370a2a7c6e4431af36aa7 Mon Sep 17 00:00:00 2001 From: Dietmar Winkler Date: Tue, 31 Mar 2026 12:04:49 +0200 Subject: [PATCH 03/10] Renamed the combined model --- .../Waterway/{RunOff_zones.mo => RunOff.mo} | 32 +++---- OpenHPL/Waterway/RunOff_zones_input.mo | 89 ------------------- OpenHPL/Waterway/package.order | 3 +- 3 files changed, 17 insertions(+), 107 deletions(-) rename OpenHPL/Waterway/{RunOff_zones.mo => RunOff.mo} (98%) delete mode 100644 OpenHPL/Waterway/RunOff_zones_input.mo diff --git a/OpenHPL/Waterway/RunOff_zones.mo b/OpenHPL/Waterway/RunOff.mo similarity index 98% rename from OpenHPL/Waterway/RunOff_zones.mo rename to OpenHPL/Waterway/RunOff.mo index 5af86e9..11d61e4 100644 --- a/OpenHPL/Waterway/RunOff_zones.mo +++ b/OpenHPL/Waterway/RunOff.mo @@ -1,21 +1,21 @@ within OpenHPL.Waterway; -model RunOff_zones "Run off model. (with 10 height zones)" +model RunOff "Run off model. (with 10 height zones)" extends OpenHPL.Icons.RunOff; - - parameter Integer N = 10 "Number of height zones" + + parameter Integer N = 10 "Number of height zones" annotation (Dialog(group = "Geometry")); - parameter SI.Area A[N] = ones(N) * 41.3e6 "Catchment area" + parameter SI.Area A[N] = ones(N) * 41.3e6 "Catchment area" annotation (Dialog(group = "Geometry")); parameter SI.PerUnit a_L[N] = {15.43, 3.97, 1.79, 0.81, 1.27, 1.44, 1.03, 2.32, 1.31, 0.57} .* 1e6 ./ A "Fractional area covered by lakes" annotation (Dialog(group = "Geometry")); - - parameter Modelica.Units.NonSI.Temperature_degC T_t=1 "Threshold temperature" + + parameter Modelica.Units.NonSI.Temperature_degC T_t=1 "Threshold temperature" annotation (Dialog(group="Physically-based parameters")); parameter Real k_m(unit="m/deg/s") = 4e-3 / 86400 "Melting factor" annotation (Dialog(group = "Physically-based parameters")); parameter SI.Volume g_T = 150e-3 "Ground saturation threshold" annotation (Dialog(group = "Physically-based parameters")); - + parameter SI.Volume s_T = 20e-3 "Soil zone saturation threshold" annotation (Dialog(group = "Empirical parameters")); parameter SI.Frequency a_1 = 0.547 / 86400 "Discharge frequency for surface runoff" @@ -34,7 +34,7 @@ model RunOff_zones "Run off model. (with 10 height zones)" annotation (Dialog(group = "Empirical parameters")); parameter Real CE(unit="1/deg") = 0.04 "Model parameter for adjusted evapotranspiration" annotation (Dialog(group = "Empirical parameters")); - + parameter Boolean useInput = false "If true, meteorological data is provided via input variables instead of data files" annotation (Dialog(tab = "Input data"), choices(checkBox = true)); parameter String fileName_temp = Modelica.Utilities.Files.loadResource("modelica://OpenHPL/Resources/Tables/Temp_var.txt") "File with temperature variations in different height zones" @@ -99,14 +99,14 @@ model RunOff_zones "Run off model. (with 10 height zones)" Modelica.Blocks.Sources.CombiTimeTable evap_var(tableOnFile = true, fileName = fileName_evap, columns = columns_evap, tableName = tableName_evap) if not useInput; Modelica.Blocks.Sources.CombiTimeTable month_temp(tableOnFile = true, fileName = fileName_month_temp, columns = columns_month_temp, tableName = tableName_month_temp) if not useInput; Modelica.Blocks.Sources.CombiTimeTable flow_var(tableOnFile = true, fileName = fileName_flow, columns = columns_flow, tableName = tableName_flow) if not useInput; - + input Real temp_input[N] if useInput "Zone temperature [degC]"; input Real prec_input[N] if useInput "Zone precipitation [mm/day]"; input Real evap_input if useInput "Potential evapotranspiration [mm/day]"; input Real month_temp_input[N] if useInput "Monthly average zone temperature [degC]"; input Real flow_input if useInput "Observed flow for R2 calculation [m3/s]"; - - Modelica.Blocks.Interfaces.RealOutput Vdot_runoff "Output connector" + + Modelica.Blocks.Interfaces.RealOutput Vdot_runoff "Output connector" annotation ( Placement(transformation(extent = {{90, -10}, {110, 10}}), iconTransformation(extent = {{80, -20}, {120, 20}}))); initial equation @@ -129,7 +129,7 @@ equation Vdot_p_r[i] = if T[i] > T_t then Vdot_p[i] * PCORR * (1 - a_L[i]) else 0; Vdot_d2w[i] = if T[i] > T_t and V_s_d[i] >= 0 then k_m * (T[i] - T_t) * (1 - a_L[i]) else 0; Vdot_s2g[i] = Vdot_p_r[i] + Vdot_d2w[i]; - + der(V_g_w[i]) = Vdot_s2g[i] - Vdot_g2s[i] - a_e[i] .* Vdot_g_e[i] "Ground zone (Soil moisure)"; Vdot_g2s[i] = if V_g_w[i] >= 0 and V_g_w[i] < g_T then (V_g_w[i] / g_T) ^ beta * Vdot_s2g[i] else Vdot_s2g[i]; if useInput then @@ -139,19 +139,19 @@ equation end if; Vdot_g_e[i] = if V_g_w[i] < g_T then V_g_w[i] / g_T * Vdot_epot[i] else Vdot_epot[i]; a_e[i] = if V_s_d[i] < err then 1 else 0; - + der(V_s_w[i]) = Vdot_g2s[i] - a_sw[i] * Vdot_s2b[i] - Vdot_s2sr[i] - Vdot_s2fr[i] "Soil zone (Upprec zone)"; Vdot_s2b[i] = (1 - a_L[i]) * PERC; Vdot_s2sr[i] = if V_s_w[i] > s_T then a_1 * (V_s_w[i] - s_T) else 0; Vdot_s2fr[i] = a_2 * V_s_w[i]; a_sw[i] = if Vdot_g2s[i] < Vdot_s2b[i] then 0 else 1; - + der(V_b_w[i]) = Vdot_s2b[i] + Vdot_pl[i] - Vdot_b2br[i] - Vdot_l_e[i] "Basement zone (Lower zone)"; Vdot_pl[i] = a_L[i] * Vdot_p[i]; Vdot_b2br[i] = a_3 * V_b_w[i]; Vdot_l_e[i] = a_L[i] * Vdot_epot[i]; end for; - + if useInput then F_o = (flow_input - 17.230144) ^ 2; F_e = (flow_input - Vdot_tot) ^ 2; @@ -371,4 +371,4 @@ recession constant for the base runoff in lower zone. Finally, the user can also the data for the CombiTimeTable models are stored.

")); -end RunOff_zones; \ No newline at end of file +end RunOff; \ No newline at end of file diff --git a/OpenHPL/Waterway/RunOff_zones_input.mo b/OpenHPL/Waterway/RunOff_zones_input.mo deleted file mode 100644 index dfd4249..0000000 --- a/OpenHPL/Waterway/RunOff_zones_input.mo +++ /dev/null @@ -1,89 +0,0 @@ -within OpenHPL.Waterway; -model RunOff_zones_input "Run off model without input data (inputs could be specified in Python)" - extends OpenHPL.Icons.RunOff; - parameter Integer N = 10 "# of height zones" annotation ( - Dialog(group = "Geometry")); - parameter Modelica.Units.NonSI.Temperature_degC T_t = 1 "Threshold temperature" annotation ( - Dialog(group = "Physically-based parameters")); - parameter SI.Area A[N] = ones(N) * 41.3e6 "Catchment area" annotation ( - Dialog(group = "Geometry")); - parameter Real s_T = 20e-3 "Soil zone saturation threshold, m" annotation ( - Dialog(group = "Empirical parameters")), a_1 = 0.547 / 86400 "Discharge frequency for surface runoff, 1/sec" annotation ( - Dialog(group = "Empirical parameters")), a_2 = 0.489 / 86400 "Discharge frequency for fast runoff, 1/sec" annotation ( - Dialog(group = "Empirical parameters")), a_3 = 0.0462 / 86400 "Discharge frequency for base runoff, 1/sec" annotation ( - Dialog(group = "Empirical parameters")), a_L[N] = {15.43, 3.97, 1.79, 0.81, 1.27, 1.44, 1.03, 2.32, 1.31, 0.57} .* 1e6 ./ A "Fractional area covered by lakes, -" annotation ( - Dialog(group = "Geometry")), g_T = 150e-3 "Ground saturation threshold, m" annotation ( - Dialog(group = "Physically-based parameters")), precC = 0.6e-3 / 86400 "preccolation from soil zone to base zone, m/sec" annotation ( - Dialog(group = "Empirical parameters")), beta = 2 "Ground zone shape coefficient, -" annotation ( - Dialog(group = "Empirical parameters")), k_m = 4e-3 / 86400 "Melting factor, m/deg/sec" annotation ( - Dialog(group = "Physically-based parameters")), PCORR = 1.05 "Precipitation correction - Rainfall, -" annotation ( - Dialog(group = "Empirical parameters")), SCORR = 1.2 "Precipitation correction - Snowfall, -" annotation ( - Dialog(group = "Empirical parameters")), CE = 0.04 "Model parameter for adjusted evapotranspiration, 1/deg" annotation ( - Dialog(group = "Empirical parameters")); - input Real temp_var[N], prec_var[N], evap_var, month_temp[N], flow_var; - SI.Height V_s_w[N] "Water content in soil zone", V_b_w[N] "Water content in base zone", V_g_w[N] "Water content in ground zone", V_s_d[N] "Dry snow"; - //V_s_s[N] "Soggy snow"; - SI.VolumeFlowRate Vdot_tot "Total runoff"; - SI.Velocity Vdot_s2b[N] "Runoff rate from soil zone to base zone", Vdot_pl[N] "Precipitation in lake", Vdot_b2br[N] "Runoff rate from base zone t obase runoff", Vdot_l_e[N] "Rate of evapotranspiration from lake", Vdot_g2s[N] "Runoff rate from ground zone to soil zone", Vdot_s2sr[N] "Runoff rate from soil zone to surface runoff", Vdot_s2fr[N] "Runoff rate from soil zone to fast runoff", Vdot_s2g[N] "Runoff rate from snow zone to ground zone", Vdot_g_e[N] "Evapotranspiration rate from ground zone", Vdot_p_r[N] "Precipitation in mainland in the form of snow", Vdot_d2w[N] "Melting rate from dry snow form to water snow form", Vdot_p_s[N] "Precipitation in mainland in the form of snow", Vdot_epot[N] "Evapotranspiration"; - //Vdot_w2d[N] "Freezing rate from water snow form to dry snow form"; - Modelica.Units.NonSI.Temperature_degC T[N] "Ambient temperature"; - SI.Velocity Vdot_p[N] "Precipitation"; - Real a_e[N], a_sw[N], F_o, F_e, R2, err = 0.5e-3 "Small error, m"; - Modelica.Blocks.Interfaces.RealOutput Vdot_runoff annotation ( - Placement(transformation(extent = {{90, -10}, {110, 10}}), iconTransformation(extent = {{80, -20}, {120, 20}}))); -initial equation - //V_s_s = zeros(N); - V_s_d = zeros(N); - V_g_w = zeros(N); - V_s_w = zeros(N); - V_b_w = zeros(N); -equation - // Total runoff - //der(Vdot_tot) = if V_s_w > s_T then (a_1*(V_s_w-s_T)+a_2*V_s_w)*(1-sum(a_L)/N)*sum(A) + sum(A)*a_3*V_b_w else a_2*V_s_w*(1-sum(a_L)/N)*sum(A) + sum(A)*a_3*V_b_w; - Vdot_tot = sum(A .* (Vdot_b2br + Vdot_s2sr + Vdot_s2fr)); - for i in 1:N loop - // Snow zone (Snow rourine) - //T[i] = temp_var.y[i]; - T[i] = temp_var[i]; - //Vdot_p[i] = prec_var.y[i] * 1e-3 / 86400; - Vdot_p[i] = prec_var[i] * 1e-3 / 86400; - der(V_s_d[i]) = Vdot_p_s[i] - Vdot_d2w[i]; - // + Vdot_w2d[i]; - //der(V_s_s[i]) = Vdot_p_r[i] - Vdot_w2d[i] + Vdot_d2w[i] - Vdot_s2g[i]; - Vdot_p_s[i] = if T[i] <= T_t then Vdot_p[i] * PCORR * SCORR * (1 - a_L[i]) else 0; - Vdot_p_r[i] = if T[i] > T_t then Vdot_p[i] * PCORR * (1 - a_L[i]) else 0; - Vdot_d2w[i] = if T[i] > T_t and V_s_d[i] >= 0 then k_m * (T[i] - T_t) * (1 - a_L[i]) else 0; - //Vdot_w2d[i] = if T[i]<=T_t and V_s_s[i]>=0 then k_m*(T_t-T[i]) else 0; - //Vdot_s2g[i] = if T[i]>T_t and V_s_s[i]=err and V_s_d[i]>=err then (1+a_w)*Vdot_d2w[i]+Vdot_p_r[i] else 0; - Vdot_s2g[i] = Vdot_p_r[i] + Vdot_d2w[i]; - // Ground zone (Soil moisure) - //der(V_g_w[i]) = Vdot_s2g[i] - Vdot_g2s[i] - (1-a)*Vdot_g_e[i]; - der(V_g_w[i]) = Vdot_s2g[i] - Vdot_g2s[i] - a_e[i] .* Vdot_g_e[i]; - Vdot_g2s[i] = if V_g_w[i] >= 0 and V_g_w[i] < g_T then (V_g_w[i] / g_T) ^ beta * Vdot_s2g[i] else Vdot_s2g[i]; - Vdot_epot[i] = evap_var * 1e-3 / 86400 * (1 + CE * (T[i] - month_temp[i])); - Vdot_g_e[i] = if V_g_w[i] < g_T then V_g_w[i] / g_T * Vdot_epot[i] else Vdot_epot[i]; - a_e[i] = if V_s_d[i] < err then 1 else 0; - // Soil zone (Upprec zone) - der(V_s_w[i]) = Vdot_g2s[i] - a_sw[i] * Vdot_s2b[i] - Vdot_s2sr[i] - Vdot_s2fr[i]; - Vdot_s2b[i] = (1 - a_L[i]) * precC; - Vdot_s2sr[i] = if V_s_w[i] > s_T then a_1 * (V_s_w[i] - s_T) else 0; - Vdot_s2fr[i] = a_2 * V_s_w[i]; - a_sw[i] = if Vdot_g2s[i] < Vdot_s2b[i] then 0 else 1; - // Basement zone (Lower zone) - der(V_b_w[i]) = Vdot_s2b[i] + Vdot_pl[i] - Vdot_b2br[i] - Vdot_l_e[i]; - Vdot_pl[i] = a_L[i] * Vdot_p[i]; - Vdot_b2br[i] = a_3 * V_b_w[i]; - Vdot_l_e[i] = a_L[i] * Vdot_epot[i]; - end for; - // Error - F_o = (flow_var - 17.230144) ^ 2; - F_e = (flow_var - Vdot_tot) ^ 2; - R2 = 1 - F_e / F_o; - Vdot_runoff = Vdot_tot; - annotation (preferredView="info", - Documentation(info=" -

This is the same hydrology model that is based on the HBV hydrological model.

-

This model can be used to define the inflow (runoff) to the reservoir. -Here the input data is not specified.

-")); -end RunOff_zones_input; diff --git a/OpenHPL/Waterway/package.order b/OpenHPL/Waterway/package.order index d380816..c39b24e 100644 --- a/OpenHPL/Waterway/package.order +++ b/OpenHPL/Waterway/package.order @@ -9,8 +9,7 @@ SurgeTank Reservoir Gate Gate_HR -RunOff_zones -RunOff_zones_input +RunOff VolumeFlowSource Penstock OpenChannel From ade7ad78d7b870f7d80667f74200709ca5b3a60e Mon Sep 17 00:00:00 2001 From: Dietmar Winkler Date: Tue, 31 Mar 2026 12:07:35 +0200 Subject: [PATCH 04/10] Fix the connector position --- OpenHPL/Waterway/RunOff.mo | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/OpenHPL/Waterway/RunOff.mo b/OpenHPL/Waterway/RunOff.mo index 11d61e4..c99b06b 100644 --- a/OpenHPL/Waterway/RunOff.mo +++ b/OpenHPL/Waterway/RunOff.mo @@ -108,7 +108,7 @@ model RunOff "Run off model. (with 10 height zones)" Modelica.Blocks.Interfaces.RealOutput Vdot_runoff "Output connector" annotation ( - Placement(transformation(extent = {{90, -10}, {110, 10}}), iconTransformation(extent = {{80, -20}, {120, 20}}))); + Placement(transformation(extent = {{100, -10}, {120, 10}}))); initial equation V_s_d = zeros(N); V_g_w = zeros(N); From 9257c3af30a414a75aa897a604c03b29bad9ffb5 Mon Sep 17 00:00:00 2001 From: Dietmar Winkler Date: Tue, 31 Mar 2026 13:43:32 +0200 Subject: [PATCH 05/10] Sync the documentation --- OpenHPL/Waterway/RunOff.mo | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/OpenHPL/Waterway/RunOff.mo b/OpenHPL/Waterway/RunOff.mo index c99b06b..e80e4fd 100644 --- a/OpenHPL/Waterway/RunOff.mo +++ b/OpenHPL/Waterway/RunOff.mo @@ -350,25 +350,34 @@ the base, quick and surface runoffs are added together.
Implementation

-This hydrology model is encoded in the OpenHPL library as the RunOff_zones unit where the main +This hydrology model is encoded in the OpenHPL library as the RunOff unit where the main defined variable is the total runoff from the catchment. This unit uses the standard Modelica connector RealOutput connector as an output from the model that can be connected to, for example, simple reservoir model Reservoir unit.

-In order to get historic information about the air temperature, precipitation, and potential evapotranspiration for each -of the elevation zones, the standard Modelica CombiTimeTable source models are used in order to read this -data from the text files. -

+Meteorological inputs (air temperature, precipitation, potential evapotranspiration, and monthly average temperature for +each elevation zone) can be provided in two ways, controlled by the boolean parameter useInput: +

+
    +
  • When useInput = false (default), historic data is read from text files using the standard Modelica +CombiTimeTable source models. The file names and table names are specified via the parameters in the +Input data tab.
    • +
  • When useInput = true, the inputs are supplied through input Real connectors +(temp_input, prec_input, evap_input, month_temp_input, and +flow_input), allowing the model to be driven by external signals at runtime.
    • +
Parameters

-When the RunOff_zones unit is in use, the user can specify the required geometry parameters for the catchment: +When the RunOff unit is in use, the user can specify the required geometry parameters for the catchment: the number of elevation zones, all hydrology parameters such as threshold temperatures, degree-day factor, precipitation correction coefficients, field capacity and β parameter in soil moisture routine, threshold level for quick runoff in upper zone, percolation from upper zone to lower zone, recession constants for the surface and quick runoffs in upper zone, and -recession constant for the base runoff in lower zone. Finally, the user can also specify the info about the text files where -the data for the CombiTimeTable models are stored. +recession constant for the base runoff in lower zone. The boolean parameter useInput selects whether +meteorological data is read from files (false) or supplied via input connectors (true). When +useInput = false, the user can also specify the text files where the data for the +CombiTimeTable models are stored.

")); end RunOff; \ No newline at end of file From ae8bcb615b8f4c7a96c64a04c20dfa247b043e0d Mon Sep 17 00:00:00 2001 From: Dietmar Winkler Date: Tue, 31 Mar 2026 13:58:07 +0200 Subject: [PATCH 06/10] Use Real interfaces instead of simple input --- OpenHPL/Waterway/RunOff.mo | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/OpenHPL/Waterway/RunOff.mo b/OpenHPL/Waterway/RunOff.mo index e80e4fd..fbebddb 100644 --- a/OpenHPL/Waterway/RunOff.mo +++ b/OpenHPL/Waterway/RunOff.mo @@ -100,11 +100,11 @@ model RunOff "Run off model. (with 10 height zones)" Modelica.Blocks.Sources.CombiTimeTable month_temp(tableOnFile = true, fileName = fileName_month_temp, columns = columns_month_temp, tableName = tableName_month_temp) if not useInput; Modelica.Blocks.Sources.CombiTimeTable flow_var(tableOnFile = true, fileName = fileName_flow, columns = columns_flow, tableName = tableName_flow) if not useInput; - input Real temp_input[N] if useInput "Zone temperature [degC]"; - input Real prec_input[N] if useInput "Zone precipitation [mm/day]"; - input Real evap_input if useInput "Potential evapotranspiration [mm/day]"; - input Real month_temp_input[N] if useInput "Monthly average zone temperature [degC]"; - input Real flow_input if useInput "Observed flow for R2 calculation [m3/s]"; + Modelica.Blocks.Interfaces.RealInput temp_input[N] if useInput "Zone temperature [degC]"; + Modelica.Blocks.Interfaces.RealInput prec_input[N] if useInput "Zone precipitation [mm/day]"; + Modelica.Blocks.Interfaces.RealInput evap_input if useInput "Potential evapotranspiration [mm/day]"; + Modelica.Blocks.Interfaces.RealInput month_temp_input[N] if useInput "Monthly average zone temperature [degC]"; + Modelica.Blocks.Interfaces.RealInput flow_input if useInput "Observed flow for R2 calculation [m3/s]"; Modelica.Blocks.Interfaces.RealOutput Vdot_runoff "Output connector" annotation ( From 2104d994495a7bf9f38beab03cb63f39c0db0c39 Mon Sep 17 00:00:00 2001 From: Dietmar Winkler Date: Tue, 31 Mar 2026 14:21:12 +0200 Subject: [PATCH 07/10] Fix balance for conditional inputs --- OpenHPL/Waterway/RunOff.mo | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/OpenHPL/Waterway/RunOff.mo b/OpenHPL/Waterway/RunOff.mo index fbebddb..cc585b5 100644 --- a/OpenHPL/Waterway/RunOff.mo +++ b/OpenHPL/Waterway/RunOff.mo @@ -94,17 +94,17 @@ model RunOff "Run off model. (with 10 height zones)" SI.VolumeFlowRate Vdot_p[N] "Precipitation"; Real a_e[N], a_sw[N], F_o, F_e, R2, err = 0.5e-3 "Small error, m"; - Modelica.Blocks.Sources.CombiTimeTable temp_var(tableOnFile = true, columns = columns_temp, tableName = tableName_temp, fileName = fileName_temp) if not useInput; - Modelica.Blocks.Sources.CombiTimeTable prec_var(tableOnFile = true, fileName = fileName_prec, columns = columns_prec, tableName = tableName_prec) if not useInput; - Modelica.Blocks.Sources.CombiTimeTable evap_var(tableOnFile = true, fileName = fileName_evap, columns = columns_evap, tableName = tableName_evap) if not useInput; - Modelica.Blocks.Sources.CombiTimeTable month_temp(tableOnFile = true, fileName = fileName_month_temp, columns = columns_month_temp, tableName = tableName_month_temp) if not useInput; - Modelica.Blocks.Sources.CombiTimeTable flow_var(tableOnFile = true, fileName = fileName_flow, columns = columns_flow, tableName = tableName_flow) if not useInput; + Modelica.Blocks.Sources.CombiTimeTable temp_var(tableOnFile = true, columns = columns_temp, tableName = tableName_temp, fileName = fileName_temp); + Modelica.Blocks.Sources.CombiTimeTable prec_var(tableOnFile = true, fileName = fileName_prec, columns = columns_prec, tableName = tableName_prec); + Modelica.Blocks.Sources.CombiTimeTable evap_var(tableOnFile = true, fileName = fileName_evap, columns = columns_evap, tableName = tableName_evap); + Modelica.Blocks.Sources.CombiTimeTable month_temp(tableOnFile = true, fileName = fileName_month_temp, columns = columns_month_temp, tableName = tableName_month_temp); + Modelica.Blocks.Sources.CombiTimeTable flow_var(tableOnFile = true, fileName = fileName_flow, columns = columns_flow, tableName = tableName_flow); - Modelica.Blocks.Interfaces.RealInput temp_input[N] if useInput "Zone temperature [degC]"; - Modelica.Blocks.Interfaces.RealInput prec_input[N] if useInput "Zone precipitation [mm/day]"; - Modelica.Blocks.Interfaces.RealInput evap_input if useInput "Potential evapotranspiration [mm/day]"; - Modelica.Blocks.Interfaces.RealInput month_temp_input[N] if useInput "Monthly average zone temperature [degC]"; - Modelica.Blocks.Interfaces.RealInput flow_input if useInput "Observed flow for R2 calculation [m3/s]"; + Modelica.Blocks.Interfaces.RealInput temp_input[N] "Zone temperature [degC]"; + Modelica.Blocks.Interfaces.RealInput prec_input[N] "Zone precipitation [mm/day]"; + Modelica.Blocks.Interfaces.RealInput evap_input "Potential evapotranspiration [mm/day]"; + Modelica.Blocks.Interfaces.RealInput month_temp_input[N] "Monthly average zone temperature [degC]"; + Modelica.Blocks.Interfaces.RealInput flow_input "Observed flow for R2 calculation [m3/s]"; Modelica.Blocks.Interfaces.RealOutput Vdot_runoff "Output connector" annotation ( @@ -159,6 +159,13 @@ equation F_o = (flow_var.y[1] - 17.230144) ^ 2; F_e = (flow_var.y[1] - Vdot_tot) ^ 2; end if; + if not useInput then + temp_input = zeros(N); + prec_input = zeros(N); + evap_input = 0; + month_temp_input = zeros(N); + flow_input = 0; + end if; R2 = 1 - F_e / F_o; Vdot_runoff = Vdot_tot; annotation (preferredView="info", From c21caedb5de7869a14e04b9059170f4469bd31a9 Mon Sep 17 00:00:00 2001 From: Dietmar Winkler Date: Tue, 31 Mar 2026 14:32:23 +0200 Subject: [PATCH 08/10] Unit fix --- OpenHPL/Waterway/RunOff.mo | 42 +++++++++++++++++++------------------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/OpenHPL/Waterway/RunOff.mo b/OpenHPL/Waterway/RunOff.mo index cc585b5..5532f8a 100644 --- a/OpenHPL/Waterway/RunOff.mo +++ b/OpenHPL/Waterway/RunOff.mo @@ -13,10 +13,10 @@ model RunOff "Run off model. (with 10 height zones)" annotation (Dialog(group="Physically-based parameters")); parameter Real k_m(unit="m/deg/s") = 4e-3 / 86400 "Melting factor" annotation (Dialog(group = "Physically-based parameters")); - parameter SI.Volume g_T = 150e-3 "Ground saturation threshold" + parameter SI.Length g_T = 150e-3 "Ground saturation threshold" annotation (Dialog(group = "Physically-based parameters")); - parameter SI.Volume s_T = 20e-3 "Soil zone saturation threshold" + parameter SI.Length s_T = 20e-3 "Soil zone saturation threshold" annotation (Dialog(group = "Empirical parameters")); parameter SI.Frequency a_1 = 0.547 / 86400 "Discharge frequency for surface runoff" annotation (Dialog(group = "Empirical parameters")); @@ -24,7 +24,7 @@ model RunOff "Run off model. (with 10 height zones)" annotation (Dialog(group = "Empirical parameters")); parameter SI.Frequency a_3 = 0.0462 / 86400 "Discharge frequency for base runoff" annotation (Dialog(group = "Empirical parameters")); - parameter SI.VolumeFlowRate PERC = 0.6e-3 / 86400 "Percolation from soil zone to base zone" + parameter SI.Velocity PERC = 0.6e-3 / 86400 "Percolation from soil zone to base zone" annotation (Dialog(group = "Empirical parameters")); parameter Real beta = 2 "Ground zone shape coefficient" annotation (Dialog(group = "Empirical parameters")); @@ -71,27 +71,27 @@ model RunOff "Run off model. (with 10 height zones)" parameter Integer columns_flow[:] = {2} "Column with real observed run off" annotation (Dialog(tab = "Input data", group = "Real run off", enable=not useInput)); - SI.Volume V_s_w[N] "Water content in soil zone"; - SI.Volume V_b_w[N] "Water content in base zone"; - SI.Volume V_g_w[N] "Water content in ground zone"; - SI.Volume V_s_d[N] "Dry snow"; + SI.Length V_s_w[N] "Water content in soil zone per Area"; + SI.Length V_b_w[N] "Water content in base zone per Area"; + SI.Length V_g_w[N] "Water content in ground zone per Area"; + SI.Length V_s_d[N] "Dry snow"; SI.VolumeFlowRate Vdot_tot "Total runoff"; - SI.VolumeFlowRate Vdot_s2b[N] "Runoff rate from soil zone to base zone"; - SI.VolumeFlowRate Vdot_pl[N] "Precipitation in lake"; - SI.VolumeFlowRate Vdot_b2br[N] "Runoff rate from base zone t obase runoff"; - SI.VolumeFlowRate Vdot_l_e[N] "Rate of evapotranspiration from lake"; - SI.VolumeFlowRate Vdot_g2s[N] "Runoff rate from ground zone to soil zone"; - SI.VolumeFlowRate Vdot_s2sr[N] "Runoff rate from soil zone to surface runoff"; - SI.VolumeFlowRate Vdot_s2fr[N] "Runoff rate from soil zone to fast runoff"; - SI.VolumeFlowRate Vdot_s2g[N] "Runoff rate from snow zone to ground zone"; - SI.VolumeFlowRate Vdot_g_e[N] "Evapotranspiration rate from ground zone"; - SI.VolumeFlowRate Vdot_p_r[N] "Precipitation in mainland in the form of snow"; - SI.VolumeFlowRate Vdot_d2w[N] "Melting rate from dry snow form to water snow form"; - SI.VolumeFlowRate Vdot_p_s[N] "Precipitation in mainland in the form of snow"; - SI.VolumeFlowRate Vdot_epot[N] "Evapotranspiration"; + SI.Velocity Vdot_s2b[N] "Runoff rate from soil zone to base zone"; + SI.Velocity Vdot_pl[N] "Precipitation in lake"; + SI.Velocity Vdot_b2br[N] "Runoff rate from base zone to base runoff"; + SI.Velocity Vdot_l_e[N] "Rate of evapotranspiration from lake"; + SI.Velocity Vdot_g2s[N] "Runoff rate from ground zone to soil zone"; + SI.Velocity Vdot_s2sr[N] "Runoff rate from soil zone to surface runoff"; + SI.Velocity Vdot_s2fr[N] "Runoff rate from soil zone to fast runoff"; + SI.Velocity Vdot_s2g[N] "Runoff rate from snow zone to ground zone"; + SI.Velocity Vdot_g_e[N] "Evapotranspiration rate from ground zone"; + SI.Velocity Vdot_p_r[N] "Precipitation in mainland in the form of snow"; + SI.Velocity Vdot_d2w[N] "Melting rate from dry snow form to water snow form"; + SI.Velocity Vdot_p_s[N] "Precipitation in mainland in the form of snow"; + SI.Velocity Vdot_epot[N] "Evapotranspiration"; Modelica.Units.NonSI.Temperature_degC T[N] "Ambient temperature"; - SI.VolumeFlowRate Vdot_p[N] "Precipitation"; + SI.Velocity Vdot_p[N] "Precipitation"; Real a_e[N], a_sw[N], F_o, F_e, R2, err = 0.5e-3 "Small error, m"; Modelica.Blocks.Sources.CombiTimeTable temp_var(tableOnFile = true, columns = columns_temp, tableName = tableName_temp, fileName = fileName_temp); From 24ee4dcb9ccd857d56ecea94eeac42a37f312c6e Mon Sep 17 00:00:00 2001 From: Dietmar Winkler Date: Tue, 31 Mar 2026 14:45:37 +0200 Subject: [PATCH 09/10] Fix SI unit --- OpenHPL/Waterway/RunOff.mo | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/OpenHPL/Waterway/RunOff.mo b/OpenHPL/Waterway/RunOff.mo index 5532f8a..d008b59 100644 --- a/OpenHPL/Waterway/RunOff.mo +++ b/OpenHPL/Waterway/RunOff.mo @@ -11,7 +11,7 @@ model RunOff "Run off model. (with 10 height zones)" parameter Modelica.Units.NonSI.Temperature_degC T_t=1 "Threshold temperature" annotation (Dialog(group="Physically-based parameters")); - parameter Real k_m(unit="m/deg/s") = 4e-3 / 86400 "Melting factor" + parameter Real k_m(unit="m/(deg.s)") = 4e-3 / 86400 "Melting factor" annotation (Dialog(group = "Physically-based parameters")); parameter SI.Length g_T = 150e-3 "Ground saturation threshold" annotation (Dialog(group = "Physically-based parameters")); @@ -387,4 +387,4 @@ meteorological data is read from files (false) or supplied via inpu CombiTimeTable models are stored.

")); -end RunOff; \ No newline at end of file +end RunOff; From f4cf1aa97968dc89ea2a278da3d8c4e15e4fabc1 Mon Sep 17 00:00:00 2001 From: Dietmar Winkler Date: Tue, 31 Mar 2026 14:59:21 +0200 Subject: [PATCH 10/10] Make use of conditional connectors for variability --- OpenHPL/Waterway/RunOff.mo | 57 ++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 30 deletions(-) diff --git a/OpenHPL/Waterway/RunOff.mo b/OpenHPL/Waterway/RunOff.mo index d008b59..05a0b7c 100644 --- a/OpenHPL/Waterway/RunOff.mo +++ b/OpenHPL/Waterway/RunOff.mo @@ -100,12 +100,18 @@ model RunOff "Run off model. (with 10 height zones)" Modelica.Blocks.Sources.CombiTimeTable month_temp(tableOnFile = true, fileName = fileName_month_temp, columns = columns_month_temp, tableName = tableName_month_temp); Modelica.Blocks.Sources.CombiTimeTable flow_var(tableOnFile = true, fileName = fileName_flow, columns = columns_flow, tableName = tableName_flow); - Modelica.Blocks.Interfaces.RealInput temp_input[N] "Zone temperature [degC]"; - Modelica.Blocks.Interfaces.RealInput prec_input[N] "Zone precipitation [mm/day]"; - Modelica.Blocks.Interfaces.RealInput evap_input "Potential evapotranspiration [mm/day]"; - Modelica.Blocks.Interfaces.RealInput month_temp_input[N] "Monthly average zone temperature [degC]"; - Modelica.Blocks.Interfaces.RealInput flow_input "Observed flow for R2 calculation [m3/s]"; - + Modelica.Blocks.Interfaces.RealInput temp_input[N] if useInput "Zone temperature [degC]"; + Modelica.Blocks.Interfaces.RealInput prec_input[N] if useInput "Zone precipitation [mm/day]"; + Modelica.Blocks.Interfaces.RealInput evap_input if useInput "Potential evapotranspiration [mm/day]"; + Modelica.Blocks.Interfaces.RealInput month_temp_input[N] if useInput "Monthly average zone temperature [degC]"; + Modelica.Blocks.Interfaces.RealInput flow_input if useInput "Observed flow for R2 calculation [m3/s]"; +protected + Modelica.Blocks.Interfaces.RealInput temp_internal[N] "Internal: zone temperature"; + Modelica.Blocks.Interfaces.RealInput prec_internal[N] "Internal: zone precipitation"; + Modelica.Blocks.Interfaces.RealInput evap_internal "Internal: potential evapotranspiration"; + Modelica.Blocks.Interfaces.RealInput month_temp_internal[N] "Internal: monthly average zone temperature"; + Modelica.Blocks.Interfaces.RealInput flow_internal "Internal: observed flow"; +public Modelica.Blocks.Interfaces.RealOutput Vdot_runoff "Output connector" annotation ( Placement(transformation(extent = {{100, -10}, {120, 10}}))); @@ -117,13 +123,8 @@ initial equation equation Vdot_tot = sum(A .* (Vdot_b2br + Vdot_s2sr + Vdot_s2fr)) "Total runoff"; for i in 1:N loop - if useInput then - T[i] = temp_input[i]; - Vdot_p[i] = prec_input[i] * 1e-3 / 86400; - else - T[i] = temp_var.y[i]; - Vdot_p[i] = prec_var.y[i] * 1e-3 / 86400; - end if; + T[i] = temp_internal[i]; + Vdot_p[i] = prec_internal[i] * 1e-3 / 86400; der(V_s_d[i]) = Vdot_p_s[i] - Vdot_d2w[i]; Vdot_p_s[i] = if T[i] <= T_t then Vdot_p[i] * PCORR * SCORR * (1 - a_L[i]) else 0; Vdot_p_r[i] = if T[i] > T_t then Vdot_p[i] * PCORR * (1 - a_L[i]) else 0; @@ -132,11 +133,7 @@ equation der(V_g_w[i]) = Vdot_s2g[i] - Vdot_g2s[i] - a_e[i] .* Vdot_g_e[i] "Ground zone (Soil moisure)"; Vdot_g2s[i] = if V_g_w[i] >= 0 and V_g_w[i] < g_T then (V_g_w[i] / g_T) ^ beta * Vdot_s2g[i] else Vdot_s2g[i]; - if useInput then - Vdot_epot[i] = evap_input * 1e-3 / 86400 * (1 + CE * (T[i] - month_temp_input[i])); - else - Vdot_epot[i] = evap_var.y[1] * 1e-3 / 86400 * (1 + CE * (T[i] - month_temp.y[i])); - end if; + Vdot_epot[i] = evap_internal * 1e-3 / 86400 * (1 + CE * (T[i] - month_temp_internal[i])); Vdot_g_e[i] = if V_g_w[i] < g_T then V_g_w[i] / g_T * Vdot_epot[i] else Vdot_epot[i]; a_e[i] = if V_s_d[i] < err then 1 else 0; @@ -152,20 +149,20 @@ equation Vdot_l_e[i] = a_L[i] * Vdot_epot[i]; end for; - if useInput then - F_o = (flow_input - 17.230144) ^ 2; - F_e = (flow_input - Vdot_tot) ^ 2; - else - F_o = (flow_var.y[1] - 17.230144) ^ 2; - F_e = (flow_var.y[1] - Vdot_tot) ^ 2; - end if; + connect(temp_input, temp_internal); + connect(prec_input, prec_internal); + connect(evap_input, evap_internal); + connect(month_temp_input, month_temp_internal); + connect(flow_input, flow_internal); if not useInput then - temp_input = zeros(N); - prec_input = zeros(N); - evap_input = 0; - month_temp_input = zeros(N); - flow_input = 0; + temp_internal = temp_var.y; + prec_internal = prec_var.y; + evap_internal = evap_var.y[1]; + month_temp_internal = month_temp.y; + flow_internal = flow_var.y[1]; end if; + F_o = (flow_internal - 17.230144) ^ 2; + F_e = (flow_internal - Vdot_tot) ^ 2; R2 = 1 - F_e / F_o; Vdot_runoff = Vdot_tot; annotation (preferredView="info",