Index: doc/wflow_sbm.rst =================================================================== diff -u -rb570d5a468e4ccd9538079eb3130b1f68700cc7f -re8ef41223499a9179f6a624980359979ba66d8bd --- doc/wflow_sbm.rst (.../wflow_sbm.rst) (revision b570d5a468e4ccd9538079eb3130b1f68700cc7f) +++ doc/wflow_sbm.rst (.../wflow_sbm.rst) (revision e8ef41223499a9179f6a624980359979ba66d8bd) @@ -144,40 +144,54 @@ Glaciers -------- -If snow modeling is enabled Glacier modelling can also be enabled by including the following three entries -in the modelparameters section: +wflow\_sbm can model glacier processes if the snow model is enabled. Glacier modelling is very close to +snow modelling and considers two main processes: glacier build-up from snow turning into firn/ice (using +the HBV-light model) and glacier melt (using a temperature degree-day model). -:: - [modelparameters] - GlacierFrac=staticmaps/GlacierFrac.map,staticmap,0.0,0 - G_TT=intbl/G_TT.tbl,tbl,0.0,1,staticmaps/GlacierFrac.map - G_Cfmax=intbl/G_Cfmax.tbl,tbl,3.0,1,staticmaps/GlacierFrac.map +The definition of glacier boundaries and initial volume is defined in two staticmaps. *GlacierFrac* is a map +that gives the fraction of each grid cell covered by a glacier as a number between zero and one. *GlacierStore* +is a state map that gives the amount of water (in mm w.e.) within the glaciers at each gridcell. Because the +glacier store (GlacierStore.map) cannot be initialized by running the model for a couple of years, a default +initial state map should be supplied by placing a GlacierStore.map file in the staticmaps directory. These two +maps are prepared from available glacier datasets. +First, a fixed fraction of the snowpack on top of the glacier is converted into ice for each timestep and added +to the glacier store using the HBV-light model (Seibert et al.,2017). This fraction, defined in the lookup table +*G_SIfrac*, typically ranges from 0.001 to 0.006. -*GlacierFrac* is a map that gives the fraction of each grid cell covered by a glacier as a number between zero and one. -Furthermore two lookup tables must be defined: *G\_TT* and *G\_Cfmax*. If the air temperature, -:math:`T_{a}`, is below :math:`G\_TT (\approx0^{o}C)` +Then, when the snowpack on top of the glacier is almost all melted (snow cover < 10 mm), glacier melt is enabled +and estimated with a degree-day model. If the air temperature, +:math:`T_{a}`, is below a certain threshold :math:`G\_TT (\approx0^{o}C)` precipitation occurs as snowfall, whereas it occurs as rainfall if :math:`T_{a}\geq G\_TT`. With this the rate of glacier melt in mm is estimated as: .. math:: - Q_{m} = cfmax(T_{a}-G\_TT)\;\;;T_{a}>G\_TT + Q_{m} = G\_Cfmax(T_{a}-G\_TT)\;\;;T_{a}>G\_TT -where :math:`Q_{m}` is the rate of glacier meltand $cfmax$ is the melting factor in :math:`mm/(^{o}C*day)`. +where :math:`Q_{m}` is the rate of glacier melt and :math:`G\_Cfmax` is the melting factor in :math:`mm/(^{o}C*day)`. +Parameters *G_TT* and *G_Cfmax* are defined in two lookup tables. *G_TT* can be taken as equal to the snow TT parameter. +Values of the melting factor normally varies from one glacier to another and some values are reported in the literature. +*G_Cfmax* can also be estimated by multiplying snow Cfmax by a factor between 1 and 2, to take into account the higher +albedo of ice compared to snow. -Accumulated snow on top of the glacier is converted to ice (and will thus become part of the glacier store) if -the total snow depth > 8300 mm. An S-curve is used to smooth the transition. A maximum of 8mm/day can be converted to ice. -.. plot:: plots/glacier-plot.py +If snow modeling is enabled, Glacier modelling can also be enabled by including the following four entries in the +modelparameters section: +:: -Becuase the glacier store (GlacierStore.map) cannot be initialized by running the model a couple of year a default -initial state map should be supplied by placing a GlacierStore.map file in the staticmaps directory. This map gives -the amount of water (in mm) within the Glaciers at each gridcell) + [modelparameters] + GlacierFrac=staticmaps/GlacierFrac.map,staticmap,0.0,0 + G_TT=intbl/G_TT.tbl,tbl,0.0,1,staticmaps/GlacierFrac.map + G_Cfmax=intbl/G_Cfmax.tbl,tbl,3.0,1,staticmaps/GlacierFrac.map + G_SIfrac=intbl/G_SIfrac.tbl,tbl,0.001,1,staticmaps/GlacierFrac.map +The initial glacier volume GlacierStore.map should also be added in the staticmaps folder. + + The rainfall interception model ------------------------------- Index: wflow/wflow_hbv.py =================================================================== diff -u -ra6b939e08b980c412cc05d2dc0bf707e7f49d2c4 -re8ef41223499a9179f6a624980359979ba66d8bd --- wflow/wflow_hbv.py (.../wflow_hbv.py) (revision a6b939e08b980c412cc05d2dc0bf707e7f49d2c4) +++ wflow/wflow_hbv.py (.../wflow_hbv.py) (revision e8ef41223499a9179f6a624980359979ba66d8bd) @@ -169,6 +169,9 @@ if hasattr(self, "ReserVoirComplexLocs"): states.append("ReservoirWaterLevel") + + if hasattr(self, "GlacierFrac"): + states.append("GlacierStore") return states @@ -337,6 +340,7 @@ pcr.setglobaloption("unittrue") self.thestep = pcr.scalar(0) + self.basetimestep = 86400 #: files to be used in case of timesries (scalar) input to the model @@ -1027,6 +1031,11 @@ self.ReservoirVolume = self.ResMaxVolume * self.ResTargetFullFrac if hasattr(self, "ReserVoirComplexLocs"): self.ReservoirWaterLevel = pcr.cover(0.0) + if hasattr(self, "GlacierFrac"): + self.GlacierStore = self.wf_readmap( + os.path.join(self.Dir, "staticmaps", "GlacierStore.map"), + 55.0 * 1000, + ) else: self.wf_resume(os.path.join(self.Dir, "instate")) @@ -1190,6 +1199,7 @@ self.SnowCover = pcr.ifthenelse(self.DrySnow > 0, pcr.scalar(1), pcr.scalar(0)) self.NrCell = pcr.areatotal(self.SnowCover, self.TopoId) + # first part of precipitation is intercepted # Interception=pcr.min(InSoil,self.ICF-self.InterceptionStorage)#: Interception in mm/timestep @@ -1219,7 +1229,32 @@ else: SnowFluxFrac = self.ZeroMap MaxFlux = self.ZeroMap + + if hasattr(self, "GlacierFrac"): + """ + Run Glacier module and add the snowpack on-top of it. + Estimate the fraction of snow turned into ice (HBV-light). + Estimate glacier melt. + glacierHBV function in wflow_lib.py + """ + self.DrySnow, self.Snow2Glacier, self.GlacierStore, self.GlacierMelt = glacierHBV( + self.GlacierFrac, + self.GlacierStore, + self.DrySnow, + self.Temperature, + self.G_TT, + self.G_Cfmax, + self.G_SIfrac, + self.timestepsecs, + self.basetimestep + ) + # Convert to mm per grid cell and add to snowmelt + self.GlacierMelt = self.GlacierMelt * self.GlacierFrac + self.FreeWater = ( + self.FreeWater + self.GlacierMelt + ) + # IntEvap=pcr.min(self.InterceptionStorage,self.PotEvaporation) #: Evaporation from interception storage # self.InterceptionStorage=self.InterceptionStorage-IntEvap Index: wflow/wflow_lib.py =================================================================== diff -u -ra9498adee6baab0a0abaa331041be8948510167b -re8ef41223499a9179f6a624980359979ba66d8bd --- wflow/wflow_lib.py (.../wflow_lib.py) (revision a9498adee6baab0a0abaa331041be8948510167b) +++ wflow/wflow_lib.py (.../wflow_lib.py) (revision e8ef41223499a9179f6a624980359979ba66d8bd) @@ -322,7 +322,62 @@ Verbose = 0 +def glacierHBV(GlacierFrac, + GlacierStore, + Snow, + Temperature, + TT, + Cfmax, + G_SIfrac, + timestepsecs, + basetimestep): + """ + Run Glacier module and add the snowpack on-top of it. + First, a fraction of the snowpack is converted into ice using the HBV-light + model (fraction between 0.001-0.005 per day). + Glacier melting is modelled using a Temperature degree factor and only + occurs if the snow cover < 10 mm. + + :ivar GlacierFrac: Fraction of wflow cell covered by glaciers + :ivar GlacierStore: Volume of the galcier in the cell in mm w.e. + :ivar Snow: Snow pack on top of Glacier + :ivar Temperature: Air temperature + :ivar TT: Temperature threshold for ice melting + :ivar Cfmax: Ice degree-day factor in mm/(°C/day) + :ivar G_SIfrac: Fraction of the snow part turned into ice each timestep + :ivar timestepsecs: Model timestep in seconds + :ivar basetimestep: Model base timestep (86 400 seconds) + + :returns: Snow,Snow2Glacier,GlacierStore,GlacierMelt, + """ + + #Fraction of the snow transformed into ice (HBV-light model) + Snow2Glacier = G_SIfrac * Snow + + Snow2Glacier = pcr.ifthenelse( + GlacierFrac > 0.0, Snow2Glacier, pcr.scalar(0.0) + ) + # Max conversion to 8mm/day + Snow2Glacier = ( + pcr.min(Snow2Glacier, 8.0) * timestepsecs / basetimestep + ) + + Snow = Snow - (Snow2Glacier * GlacierFrac) + GlacierStore = GlacierStore + Snow2Glacier + + PotMelt = pcr.ifthenelse( + Temperature > TT, Cfmax * (Temperature - TT), pcr.scalar(0.0) + ) # Potential snow melt, based on temperature + + GlacierMelt = pcr.ifthenelse( + Snow < 10.0, pcr.min(PotMelt, GlacierStore), pcr.cover(0.0) + ) # actual Glacier melt + GlacierStore = GlacierStore - GlacierMelt # dry snow content + + return Snow, Snow2Glacier, GlacierStore, GlacierMelt + + def lddcreate_save( lddname, dem, Index: wflow/wflow_sbm.py =================================================================== diff -u -r11471548f856592485a77ca37df4ad838eed2eaf -re8ef41223499a9179f6a624980359979ba66d8bd --- wflow/wflow_sbm.py (.../wflow_sbm.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) +++ wflow/wflow_sbm.py (.../wflow_sbm.py) (revision e8ef41223499a9179f6a624980359979ba66d8bd) @@ -383,31 +383,6 @@ return Snow, SnowWater, SnowMelt, RainFall, SnowFall -def GlacierMelt(GlacierStore, Snow, Temperature, TT, Cfmax): - """ - Glacier modelling using a Temperature degree factor. Melting - only occurs if the snow cover > 10 mm - - - :ivar GlacierStore: - :ivar Snow: Snow pack on top of Glacier - :ivar Temperature: - - :returns: GlacierStore,GlacierMelt, - """ - - PotMelt = pcr.ifthenelse( - Temperature > TT, Cfmax * (Temperature - TT), pcr.scalar(0.0) - ) # Potential snow melt, based on temperature - - GlacierMelt = pcr.ifthenelse( - Snow > 10.0, pcr.min(PotMelt, GlacierStore), pcr.cover(0.0) - ) # actual Glacier melt - GlacierStore = GlacierStore - GlacierMelt # dry snow content - - return GlacierStore, GlacierMelt - - class WflowModel(pcraster.framework.DynamicModel): """ .. versionchanged:: 0.91 @@ -1886,33 +1861,21 @@ if hasattr(self, "GlacierFrac"): """ Run Glacier module and add the snowpack on-top of it. - Snow becomes ice when pressure is about 830 k/m^2, e.g 8300 mm - If below that a max amount of 2mm/day can be converted to glacier-ice + Estimate the fraction of snow turned into ice (HBV-light). + Estimate glacier melt. + glacierHBV function in wflow_lib.py """ - # TODO: document glacier module - self.snowdist = sCurve(self.Snow, a=8300.0, c=0.06) - self.Snow2Glacier = pcr.ifthenelse( - self.Snow > 8300, self.snowdist * (self.Snow - 8300), self.ZeroMap - ) - self.Snow2Glacier = pcr.ifthenelse( - self.GlacierFrac > 0.0, self.Snow2Glacier, self.ZeroMap - ) - # Max conversion to 8mm/day - self.Snow2Glacier = ( - pcr.min(self.Snow2Glacier, 8.0) - * self.timestepsecs - / self.basetimestep - ) - - self.Snow = self.Snow - (self.Snow2Glacier * self.GlacierFrac) - - self.GlacierStore, self.GlacierMelt = GlacierMelt( - self.GlacierStore + self.Snow2Glacier, + self.Snow, self.Snow2Glacier, self.GlacierStore, self.GlacierMelt = glacierHBV( + self.GlacierFrac, + self.GlacierStore, self.Snow, self.Temperature, self.G_TT, self.G_Cfmax, + self.G_SIfrac, + self.timestepsecs, + self.basetimestep ) # Convert to mm per grid cell and add to snowmelt self.GlacierMelt = self.GlacierMelt * self.GlacierFrac