updating snow calculation in HBV

Author: anatappe

MATILDA/MATILDA_slim/MATILDA.py

@@ -509,7 +509,7 @@ def glacier_change(output_DDM, lookup_table, glacier_profile, parameter):
     in the unit mm / day.
     """
 
-    def hbv_simulation(input_df_catchment, parameter):
+    def hbv_simulation(input_df_catchment, parameter, glacier_area=None):
         print("Running HBV routine")
         # 1. new temporary dataframe from input with daily values
         if "PE" in input_df_catchment.columns:
@@ -573,7 +573,21 @@ def hbv_simulation(input_df_catchment, parameter):
         SNOW_cal = snowfrac_cal * Prec_cal
         # snow correction factor
         SNOW_cal = parameter.SFCF * SNOW_cal
-        SNOW_cal =SNOW_cal * (1-(parameter.area_glac / parameter.area_cat))
+        # get the new glacier area for each year
+        if glacier_area is not None:
+            glacier_area = glacier_area.iloc[1:, :]
+            glacier_area["time"] = glacier_area["time"].astype(str).astype(float).astype(int)
+            SNOW2 = pd.DataFrame(SNOW_cal)
+            SNOW2["area"] = 0
+            for year in range(len(glacier_area)):
+                SNOW2["area"] = np.where(SNOW2.index.year == glacier_area["time"].iloc[year],
+                                         glacier_area["glacier_area"].iloc[year],
+                                         SNOW2["area"])
+
+            SNOW2["T2"] = SNOW2["T2"] * (1 - (SNOW2["area"] / parameter.area_cat))
+            SNOW_cal = SNOW2['T2'].squeeze()
+        else:
+            SNOW_cal["T2"] = SNOW_cal["T2"] * (1 - (parameter.area_glac / parameter.area_cat))
         # evaporation correction
         # a. calculate long-term averages of daily temperature
         Temp_mean_cal = np.array([Temp_cal.loc[Temp_cal.index.dayofyear == x].mean() \
@@ -676,8 +690,21 @@ def hbv_simulation(input_df_catchment, parameter):
         SNOW = parameter.SFCF * SNOW
         # snow correction factor
         SNOW = parameter.SFCF * SNOW
-        SNOW = SNOW * (1-(parameter.area_glac / parameter.area_cat))
-        # evaporation correction
+        # get the new glacier area for each year
+        if glacier_area is not None:
+            glacier_area = glacier_area.iloc[1:, :]
+            glacier_area["time"] = glacier_area["time"].astype(str).astype(float).astype(int)
+            SNOW2 = pd.DataFrame(SNOW)
+            SNOW2["area"] = 0
+            for year in range(len(glacier_area)):
+                SNOW2["area"] = np.where(SNOW2.index.year == glacier_area["time"].iloc[year],
+                                         glacier_area["glacier_area"].iloc[year],
+                                         SNOW2["area"])
+
+            SNOW2["T2"] = SNOW2["T2"] * (1 - (SNOW2["area"] / parameter.area_cat))
+            SNOW = SNOW2['T2'].squeeze()
+        else:
+            SNOW["T2"] = SNOW["T2"] * (1 - (parameter.area_glac / parameter.area_cat))        # evaporation correction
         # a. calculate long-term averages of daily temperature
         Temp_mean = np.array([Temp.loc[Temp.index.dayofyear == x].mean() \
                               for x in range(1, 367)])
@@ -815,7 +842,10 @@ def hbv_simulation(input_df_catchment, parameter):
         print("Finished HBV routine")
         return hbv_results
 
-    output_HBV = hbv_simulation(input_df_catchment, parameter)
+    if glacier_profile is not None:
+        output_HBV = hbv_simulation(input_df_catchment, parameter, glacier_area=glacier_change_area)
+    else:
+        output_HBV = hbv_simulation(input_df_catchment, parameter)
     output_HBV = output_HBV[parameter.sim_start:parameter.sim_end]
 
     # Output postprocessing

MATILDA/run_MATILDA.py

@@ -54,17 +54,4 @@
 output_MATILDA[4].show()
 
 output_MATILDA[0].Q_Total
-<<<<<<< HEAD
-=======
 
-##
-obs = pd.read_csv("/home/ana/Seafile/Ana-Lena_Phillip/data/input_output/input/observations/bash_kaindy/runoff_bashkaindy_04_2019-11_2020_temp_limit.csv")
-obs_poly = pd.read_csv("/home/ana/Seafile/Masterarbeit/Data/Input/discharge_bashkaingdy_polyfitted_2019-04_11-2020.csv")
-
-plt.plot(obs["Qobs"], label="Obs old (mm?)")
-#plt.plot(obs_poly["discharge_poly"], label="Obs poly (m3/s)")
-#plt.plot((obs["Qobs"] / 86400*(46.232*1000000)/1000), label="Obs (old) calc from mm to m/3")
-plt.plot((obs_poly["discharge_poly"] * 86400 / (46.232 * 1000000) * 1000), label="Obs poly from m/3 to mm")
-plt.legend()
-plt.show()
->>>>>>> 235dfbba2e4a1bcf29f4d78aee722abb240fd595

Test_area/MATILDA_hourly.py

@@ -3,7 +3,7 @@
 MATILDA (Modeling wATer resources In gLacierizeD cAtchments) is a combination of a degree day model and the HBV model (Bergstöm 1976) to compute total runoff of glacierized catchments.
 This file may use the input files created by the COSIPY-utility "aws2cosipy" as forcing data and or a simple dataframe with temperature, precipitation and if possible evapotranspiration and additional observation runoff data to validate it.
 """
-# Import all necessary python packages
+## Import all necessary python packages
 import xarray as xr
 import numpy as np
 import pandas as pd
@@ -16,7 +16,7 @@
 from matplotlib.offsetbox import AnchoredText
 from MATILDA_slim import MATILDA
 
-
+##
 # Setting the parameter for the MATILDA simulation
 def MATILDA_parameter(input_df, set_up_start=None, set_up_end=None, sim_start=None, sim_end=None, freq="D",
                       lat= None, area_cat=None, area_glac=None, ele_dat=None, ele_glac=None, ele_cat=None, parameter_df = None,
@@ -1152,9 +1152,10 @@ def MATILDA_simulation(input_df, obs=None, glacier_profile=None, output=None, se
 
 parameter = MATILDA.MATILDA_parameter(df, set_up_start='2018-01-01 00:00:00', set_up_end='2018-12-31 23:00:00',
                                       sim_start='2019-01-01 00:00:00', sim_end='2020-11-01 23:00:00', freq="D",
-                                      lat=41, area_cat=46.23, area_glac=2.566, ele_dat=2250, ele_glac=4035,
+                                      lat=41, area_cat=46.23, area_glac=2.566, ele_dat=3864, ele_glac=4035,
                                       ele_cat=3485, CFMAX_ice=5.5, CFMAX_snow=3)
 
+
 # load the right package for the right resolution
 input_df_glacier = df[parameter.sim_start:parameter.sim_end]
 if parameter.area_glac > 0: