Source code for cwatm.hydrological_modules.evaporationPot

# -------------------------------------------------------------------------
# Name:        POTENTIAL REFERENCE EVAPO(TRANSPI)RATION
# Purpose:
#
# Author:      PB
#
# Created:     10/01/2017
# Copyright:   (c) PB 2017
# -------------------------------------------------------------------------

from cwatm.management_modules.data_handling import *


[docs]class evaporationPot(object): """ POTENTIAL REFERENCE EVAPO(TRANSPI)RATION Calculate potential evapotranspiration from climate data mainly based on FAO 56 and LISVAP Based on Penman Monteith References: http://www.fao.org/docrep/X0490E/x0490e08.htm#penman%20monteith%20equation http://www.fao.org/docrep/X0490E/x0490e06.htm http://www.fao.org/docrep/X0490E/x0490e06.htm https://ec.europa.eu/jrc/en/publication/eur-scientific-and-technical-research-reports/lisvap-evaporation-pre-processor-lisflood-water-balance-and-flood-simulation-model **Global variables** ===================================== ====================================================================== ===== Variable [self.var] Description Unit ===================================== ====================================================================== ===== cropCorrect calibration factor of crop KC factor -- pet_modus Flag: index which ETP approach is used e.g. 1 for Penman-Monteith bool AlbedoCanopy Albedo of vegetation canopy (FAO,1998) default =0.23 -- AlbedoSoil Albedo of bare soil surface (Supit et. al. 1994) default = 0.15 -- AlbedoWater Albedo of water surface (Supit et. al. 1994) default = 0.05 -- dem -- lat -- co2 -- albedoLand albedo from land surface (from GlobAlbedo database) -- albedoOpenWater albedo from open water surface (from GlobAlbedo database) -- ETRef potential evapotranspiration rate from reference crop m only_radition -- TMin minimum air temperature K TMax maximum air temperature K Tavg Input, average air Temperature K Rsds short wave downward surface radiation fluxes W/m2 EAct -- Psurf Instantaneous surface pressure Pa Qair specific humidity kg/kg Rsdl long wave downward surface radiation fluxes W/m2 Wind wind speed m/s EWRef potential evaporation rate from water surface m ===================================== ====================================================================== ===== **Functions** """ def __init__(self, model): """ The constructor evaporationPot """ self.var = model.var self.model = model
[docs] def initial(self): """ Initial part of evaporation type module Load inictial parameters Note: Only run if *calc_evaporation* is True """ #self.var.sumETRef = globals.inZero.copy() self.var.cropCorrect = loadmap('crop_correct') if checkOption('calc_evaporation'): # Default calculation method is Penman Monteith # if PET_modus is missing use Penman Monteith self.var.pet_modus = 1 if "PET_modus" in option: self.var.pet_modus = checkOption('PET_modus') self.initial_1()
#if self.var.pet_modus == 1: self.initial_1() # --------------------------------------------------------------------------
[docs] def initial_1(self): """ Initial part of evaporation type module Load initial parameters 1: Penman Monteith 2: Milly and Dunne method P. C. D. Milly* and K. A. Dunne, 2016: Potential evapotranspiration and continental drying, Nature Climate Change, DOI: 10.1038/NCLIMATE3046 Energy only PET: ET=0.8(Rn ?) equation 8 3: Yang et al. Penman Montheith correction method Yang, Y., Roderick, M. L., Zhang, S., McVicar, T. R., and Donohue, R. J.: Hydrologic implications of vegetation response to elevated CO2 in climate projections, Nat. Clim. Change, 9, 44-48, 10.1038/s41558-018-0361-0, 2019. Equation 14: where the term 2.14 accounts for changing [CO2] on rs """ self.var.AlbedoCanopy = loadmap('AlbedoCanopy') self.var.AlbedoSoil = loadmap('AlbedoSoil') self.var.AlbedoWater = loadmap('AlbedoWater') if self.var.only_radition: self.var.dem = loadmap('dem') self.var.lat = loadmap('latitude')
# -------------------------------------------------------------------------- # --------------------------------------------------------------------------
[docs] def dynamic(self): """ Dynamic part of the potential evaporation module :return: - ETRef - potential reference evapotranspiration rate [m/day] - EWRef - potential evaporation rate from water surface [m/day] """ if checkOption('calc_evaporation'): if self.var.pet_modus == 1: self.dynamic_1() if self.var.pet_modus == 2: self.dynamic_2() if self.var.pet_modus == 3: self.dynamic_1() if self.var.pet_modus == 4: self.dynamic_4()
# --------------------------------------------------------------------------
[docs] def dynamic_1(self): """ Dynamic part of the potential evaporation module Based on Penman Monteith - FAO 56 """ if self.var.pet_modus == 3: #Yang et al.Penman Montheith correction method # loading CO2 concentration RCP2.6-RCP8.5 if dateVar['newStart'] or dateVar['newYear']: self.var.co2 = readnetcdf2('co2conc', dateVar['currDate'], "yearly", value="CO2", cut = False, compress= False) ESatmin = 0.6108* np.exp((17.27 * self.var.TMin) / (self.var.TMin + 237.3)) ESatmax = 0.6108* np.exp((17.27 * self.var.TMax) / (self.var.TMax + 237.3)) ESat = (ESatmin + ESatmax) / 2.0 # [KPa] # http://www.fao.org/docrep/X0490E/x0490e07.htm equation 11/12 RNup = 4.903E-9 * (((self.var.TMin + 273.16) ** 4) + ((self.var.TMax + 273.16) ** 4)) / 2 # Up longwave radiation [MJ/m2/day] LatHeatVap = 2.501 - 0.002361 * self.var.Tavg # latent heat of vaporization [MJ/kg] # -------------------------------- # if only daily calculate radiation is given instead of longwave down and shortwave down radiation if self.var.only_radition: # FAO 56 - https://www.fao.org/3/x0490E/x0490e07.htm#solar%20radiation equation 39 radian = np.pi / 180 * self.var.lat distanceSun = 1 + 0.033 * np.cos(2 * np.pi * dateVar['doy'] / 365) # Chapter 3: equation 24 declin = 0.409 * np.sin(2 * np.pi * dateVar['doy'] / 365 - 1.39) ws = np.arccos(-np.tan(radian * np.tan(declin))) Ra = 24 *60 / np.pi * 0.082 * distanceSun * (ws * np.sin(radian) * np.sin(declin) + np.cos(radian) * np.cos(declin) * np.sin(ws)) # Equation 21 Chapter 3 Rso = Ra * (0.75 + (2 * 10 ** -5 * self.var.dem)) # in MJ/m2/day # Equation 37 Chapter 3 RsRso = 1.35 * self.var.Rsds/Rso - 0.35 RsRso = np.minimum(np.maximum(RsRso, 0.05), 1) EmNet = (0.34 - 0.14 * np.sqrt(self.var.EAct)) # Eact in hPa but needed in kPa : kpa = 0.1 * hPa - conversion done in readmeteo RLN = RNup * EmNet * RsRso # Equation 39 Chapter 3 Psycon = 0.00163 * (101.3 / LatHeatVap) # psychrometric constant at sea level [mbar/deg C] #Psycon = 0.665E-3 * self.var.Psurf # psychrometric constant [kPa C-1] # http://www.fao.org/docrep/X0490E/x0490e07.htm Equation 8 # see http://www.fao.org/docrep/X0490E/x0490e08.htm#penman%20monteith%20equation Psycon = Psycon * ((293 - 0.0065 * self.var.dem) / 293) ** 5.26 # in [KPa deg C-1] # http://www.fao.org/docrep/X0490E/x0490e07.htm Equation 7 else: Psycon = 0.665E-3 * self.var.Psurf # psychrometric constant [kPa C-1] # http://www.fao.org/docrep/X0490E/x0490e07.htm Equation 8 # see http://www.fao.org/docrep/X0490E/x0490e08.htm#penman%20monteith%20equation # calculate vapor pressure # Fao 56 Page 36 # calculate actual vapour pressure if returnBool('useHuss'): # if specific humidity calculate actual vapour pressure this way self.var.EAct = (self.var.Psurf * self.var.Qair) / ((0.378 * self.var.Qair) + 0.622) # http://www.eol.ucar.edu/projects/ceop/dm/documents/refdata_report/eqns.html # (self.var.Psurf * self.var.Qair)/0.622 # old calculation not completely ok else: # if relative humidity self.var.EAct = ESat * self.var.Qair / 100.0 # longwave radiation balance RLN = RNup - self.var.Rsdl # RDL is stored on disk as W/m2 but converted in MJ/m2/s in readmeteo.py # ************************************************************ # ***** NET ABSORBED RADIATION ******************************* # ************************************************************ if returnBool('albedo'): if dateVar['newStart'] or dateVar['newMonth']: # loading every month a new map self.var.albedoLand = readnetcdf2('albedoMaps', dateVar['currDate'], useDaily='month',value='albedoLand') self.var.albedoOpenWater = readnetcdf2('albedoMaps', dateVar['currDate'], useDaily='month',value='albedoWater') RNA = np.maximum(((1 - self.var.albedoLand) * self.var.Rsds - RLN) / LatHeatVap, 0.0) RNAWater = np.maximum(((1 - self.var.albedoOpenWater) * self.var.Rsds - RLN) / LatHeatVap, 0.0) else: RNA = np.maximum(((1 - self.var.AlbedoCanopy) * self.var.Rsds - RLN) / LatHeatVap, 0.0) # net absorbed radiation of reference vegetation canopy [mm/d] # RNASoil = np.maximum(((1 - self.var.AlbedoSoil) * self.var.Rsds - RLN) / LatHeatVap, 0.0) # net absorbed radiation of bare soil surface RNAWater = np.maximum(((1 - self.var.AlbedoWater) * self.var.Rsds - RLN) / LatHeatVap, 0.0) # net absorbed radiation of water surface VapPressDef = np.maximum(ESat - self.var.EAct, 0.0) Delta = ((4098.0 * ESat) / ((self.var.Tavg + 237.3)**2)) # slope of saturated vapour pressure curve [kPa/deg C] # Equation 13 Chapter 3 # Chapter 2 Equation 6 windpart = 900 * self.var.Wind / (self.var.Tavg + 273.16) if self.var.pet_modus == 1: denominator = Delta + Psycon *(1 + 0.34 * self.var.Wind) else: # Yang et al.Penman Montheith correction method: term 2 accounts for changing [CO2] on rs. denominator = Delta + Psycon * (1 + self.var.Wind*(0.34+0.00024*(self.var.co2-300.))) numerator1 = Delta / denominator numerator2 = Psycon / denominator # the 0.408 constant is replace by 1/LatHeatVap see above RNAN = RNA * numerator1 #RNANSoil = RNASoil * numerator1 RNANWater = RNAWater * numerator1 EA = windpart * VapPressDef * numerator2 # Potential evapo(transpi)ration is calculated for two reference surfaces: # 1. Reference vegetation canopy # 2. Open water surface self.var.ETRef = (RNAN + EA) * 0.001 # potential reference evapotranspiration rate [m/day] # from mm to m with 0.001 #self.var.ESRef = RNANSoil + EA # potential evaporation rate from a bare soil surface [m/day] self.var.EWRef = (RNANWater + EA) * 0.001 ii =1
# potential evaporation rate from water surface [m/day] # -> here we are at ET0 (see http://www.fao.org/docrep/X0490E/x0490e04.htm#TopOfPage figure 4:) #self.var.sumETRef = self.var.sumETRef + self.var.ETRef*1000 #if dateVar['curr'] ==32: #ii=1 #report(decompress(self.var.sumETRef), "C:\work\output2/sumetref.map")
[docs] def dynamic_2(self): """ Dynamic part of the potential evaporation module 2: Milly and Dunne method P. C. D. Milly* and K. A. Dunne, 2016: Potential evapotranspiration and continental drying, Nature Climate Change, DOI: 10.1038/NCLIMATE3046 Energy only PET = 0.8(Rn ? ) equation 8 """ LatHeatVap = 2.501 - 0.002361 * self.var.Tavg # latent heat of vaporization [MJ/kg] RNup = 4.903E-9 * (((self.var.TMin + 273.16) ** 4) + ((self.var.TMax + 273.16) ** 4)) / 2 # Up longwave radiation [MJ/m2/day] if self.var.only_radition: # FAO 56 - https://www.fao.org/3/x0490E/x0490e07.htm#solar%20radiation equation 39 a = dateVar['doy'] #radian = np.pi / 180 * self.var.lat radian = np.pi / 180 * -20 #distanceSun = 1 + 0.033 * np.cos(2 * np.pi * dateVar['doy'] / 365) distanceSun = 1 + 0.033 * np.cos(2 * np.pi * 246 / 365) #declin = 0.409 * np.sin(2 * np.pi * dateVar['doy'] / 365 - 1.39) declin = 0.409 * np.sin(2 * np.pi * 246 / 365 - 1.39) ws = np.arccos(-np.tan(radian * np.tan(declin))) Ra = 24 *60 / np.pi * 0.082 * distanceSun * (ws * np.sin(radian) * np.sin(declin) + np.cos(radian) * np.cos(declin) * np.sin(ws)) Rso = Ra * (0.75 + (2 * 10 ** -5 * self.var.dem)) # in MJ/m2/day RsRso = 1.35 * self.var.Rsds/Rso - 0.35 RsRso = np.minimum(np.maximum(RsRso, 0.05), 1) EmNet = (0.34 - 0.14 * np.sqrt(self.var.EAct)) # Eact in hPa but needed in kPa : kpa = 0.1 * hPa - conversion done in readmeteo RLN = RNup * EmNet * RsRso else: RLN = RNup - self.var.Rsdl # RDL is stored on disk as W/m2 but converted in MJ/m2/s in readmeteo.py if returnBool('albedo'): if dateVar['newStart'] or dateVar['newMonth']: # loading every month a new map self.var.albedoLand = readnetcdf2('albedoMaps', dateVar['currDate'], useDaily='month',value='albedoLand') self.var.albedoOpenWater = readnetcdf2('albedoMaps', dateVar['currDate'], useDaily='month',value='albedoWater') RNA = np.maximum(((1 - self.var.albedoLand) * self.var.Rsds - RLN) / LatHeatVap, 0.0) RNAWater = np.maximum(((1 - self.var.albedoOpenWater) * self.var.Rsds - RLN) / LatHeatVap, 0.0) else: RNA = np.maximum(((1 - self.var.AlbedoCanopy) * self.var.Rsds - RLN) / LatHeatVap, 0.0) # net absorbed radiation of reference vegetation canopy [mm/d] # RNASoil = np.maximum(((1 - self.var.AlbedoSoil) * self.var.Rsds - RLN) / LatHeatVap, 0.0) # net absorbed radiation of bare soil surface RNAWater = np.maximum(((1 - self.var.AlbedoWater) * self.var.Rsds - RLN) / LatHeatVap, 0.0) # net absorbed radiation of water surface # 1. Reference vegetation canopy # 2. Open water surface self.var.ETRef = 0.8 * RNA * 0.001 # potential reference evapotranspiration rate [m/day] # from mm to m with 0.001 # potential evaporation rate from a bare soil surface [m/day] self.var.EWRef = 0.8 * RNAWater * 0.001
[docs] def dynamic_4(self): """ Dynamic part of the potential evaporation module 4. Priestley-Taylor 1.26 * delat https://wetlandscapes.github.io/blog/blog/penman-monteith-and-priestley-taylor/ uses only tmin, tmax, tavg, rsds, rlds (or rsd) """ ESatmin = 0.6108* np.exp((17.27 * self.var.TMin) / (self.var.TMin + 237.3)) ESatmax = 0.6108* np.exp((17.27 * self.var.TMax) / (self.var.TMax + 237.3)) ESat = (ESatmin + ESatmax) / 2.0 # [KPa] # http://www.fao.org/docrep/X0490E/x0490e07.htm equation 11/12 RNup = 4.903E-9 * (((self.var.TMin + 273.16) ** 4) + ((self.var.TMax + 273.16) ** 4)) / 2 # Up longwave radiation [MJ/m2/day] LatHeatVap = 2.501 - 0.002361 * self.var.Tavg # latent heat of vaporization [MJ/kg] # if only daily calculate radiation is given instead of longwave down and shortwave down radiation if self.var.only_radition: # FAO 56 - https://www.fao.org/3/x0490E/x0490e07.htm#solar%20radiation equation 39 a = dateVar['doy'] #radian = np.pi / 180 * self.var.lat radian = np.pi / 180 * -20 #distanceSun = 1 + 0.033 * np.cos(2 * np.pi * dateVar['doy'] / 365) distanceSun = 1 + 0.033 * np.cos(2 * np.pi * 246 / 365) #declin = 0.409 * np.sin(2 * np.pi * dateVar['doy'] / 365 - 1.39) declin = 0.409 * np.sin(2 * np.pi * 246 / 365 - 1.39) ws = np.arccos(-np.tan(radian * np.tan(declin))) Ra = 24 *60 / np.pi * 0.082 * distanceSun * (ws * np.sin(radian) * np.sin(declin) + np.cos(radian) * np.cos(declin) * np.sin(ws)) Rso = Ra * (0.75 + (2 * 10 ** -5 * self.var.dem)) # in MJ/m2/day RsRso = 1.35 * self.var.Rsds/Rso - 0.35 RsRso = np.minimum(np.maximum(RsRso, 0.05), 1) EmNet = (0.34 - 0.14 * np.sqrt(self.var.EAct)) # Eact in hPa but needed in kPa : kpa = 0.1 * hPa - conversion done in readmeteo RLN = RNup * EmNet * RsRso Psycon = 0.00163 * (101.3 / LatHeatVap) # psychrometric constant at sea level [mbar/deg C] #Psycon = 0.665E-3 * self.var.Psurf # psychrometric constant [kPa C-1] # http://www.fao.org/docrep/X0490E/x0490e07.htm Equation 8 # see http://www.fao.org/docrep/X0490E/x0490e08.htm#penman%20monteith%20equation Psycon = Psycon * ((293 - 0.0065 * self.var.dem) / 293) ** 5.26 # in [KPa deg C-1] # http://www.fao.org/docrep/X0490E/x0490e07.htm Equation 7 else: RLN = RNup - self.var.Rsdl Psycon = 0.665E-3 * self.var.Psurf # psychrometric constant [kPa C-1] # http://www.fao.org/docrep/ X0490E/ x0490e07.htm Equation 8 # see http://www.fao.org/docrep/X0490E/x0490e08.htm#penman%20monteith%20equation # RDL is stored on disk as W/m2 but converted in MJ/m2/s in readmeteo.py RNA = np.maximum(((1 - self.var.AlbedoCanopy) * self.var.Rsds - RLN) / LatHeatVap, 0.0) RNAWater = np.maximum(((1 - self.var.AlbedoWater) * self.var.Rsds - RLN) / LatHeatVap, 0.0) Delta = ((4098.0 * ESat) / ((self.var.Tavg + 237.3)**2)) # slope of saturated vapour pressure curve [mbar/deg C] RNAN = 1.1 * Delta / (Delta + Psycon) * RNA #RNANSoil = RNASoil * numerator1 RNANWater = 1.26 * Delta / (Delta + Psycon) * RNAWater # 1. Reference vegetation canopy # 2. Open water surface self.var.ETRef = RNAN * 0.001 # potential reference evapotranspiration rate [m/day] # from mm to m with 0.001 # potential evaporation rate from a bare soil surface [m/day] self.var.EWRef = RNANWater * 0.001