""" Created on March 14, 2017 Originally written by Scott Havens in 2015 @author: Micah Johnson **Creating Custom NASDE Models** -------------------------------- When creating a new NASDE model make sure you adhere to the following: 1. Add a new method with the other models with a unique name ideally with some reference to the origin of the model. For example see :func:`~smrf.envphys.snow.susong1999`. #. Add the new model to the dictionary :data:`~smrf.envphys.snow.available_models` at the bottom of this module so that :func:`~smrf.envphys.snow.calc_phase_and_density` can see it. #. Create a custom distribution function with a unique in :func:`~smrf.distribute.precipitation.distribute` to create the structure for the new model. For an example see :func:`~smrf.distribute.precipitation.distribute_for_susong1999`. #. Update documentation and run smrf! """ import numpy as np import pandas as pd [docs]def calc_phase_and_density(temperature, precipitation, nasde_model): """ Uses various new accumulated snow density models to estimate the snow density of precipitation that falls during sub-zero conditions. The models all are based on the dew point temperature and the amount of precipitation, All models used here must return a dictionary containing the keywords pcs and rho_s for percent snow and snow density respectively. Args: temperature: a single timestep of the distributed dew point temperature precipitation: a numpy array of the distributed precipitation nasde_model: string value set in the configuration file representing the method for estimating density of new snow that has just fallen. Returns: tuple: Returns a tuple containing the snow density field and the percent snow as determined by the NASDE model. - **snow_density** (*numpy.array*) - Snow density values in kg/m^3 - **perc_snow** (*numpy.array*) - Percent of the precip that is snow in values 0.0-1.0. """ # convert the inputs to numpy arrays precip = np.array(precipitation) temperature = np.array(temperature) snow_density = np.zeros(precip.shape) perc_snow = np.zeros(precip.shape) #New accumulated snow point models can go here. if nasde_model in available_models.keys(): result = available_models[nasde_model](temperature,precip) else: raise ValueError("{0} is not an implemented new accumulated snow density (NASDE) model! Check the config file under precipitation".format(nasde_model)) snow_density = result['rho_s'] perc_snow = result['pcs'] return snow_density,perc_snow [docs]def calc_perc_snow(Tpp, Tmax=0.0, Tmin=-10.0): """ Calculates the percent snow for the nasde_models piecewise_susong1999 and marks2017. Args: Tpp: A numpy array of temperature, use dew point temperature if available [degree C]. Tmax: Max temperature that the percent snow is estimated. Default is 0.0 Degrees C. Tmin: Minimum temperature that percent snow is changed. Default is -10.0 Degrees C. Returns: numpy.array: A fraction of the precip at each pixel that is snow provided by Tpp. """ #Coefficients for snow relationship Tr0 = 0.5 Pcr0 = 0.25 Pc0 = 0.75 pcs = np.zeros(Tpp.shape) pcs[Tpp <= -0.5] = 1.0 ind = (Tpp > -0.5) & (Tpp <= 0.0) if np.any(ind): pcs[ind] = (-Tpp[ind] / Tr0) * Pcr0 + Pc0 ind = (Tpp > 0.0) & (Tpp <= Tmax +1.0) if np.any(ind): pcs[ind] = (-Tpp[ind] / (Tmax + 1.0)) * Pc0 + Pc0 return pcs [docs]def check_temperature(Tpp, Tmax = 0.0, Tmin = -10.0): """ Sets the precipitation temperature and snow temperature. Args: Tpp: A numpy array of temperature, use dew point temperature if available [degrees C]. Tmax: Thresholds the max temperature of the snow [degrees C]. Tmin: Minimum temperature that the precipitation temperature [degrees C]. Returns: tuple: - **Tpp** (*numpy.array*) - Modified precipitation temperature that is thresholded with a minimum set by tmin. - **tsnow** (*numpy.array*) - Temperature of the surface of the snow set by the precipitation temperature and thresholded by tmax where tsnow > tmax = tmax. """ Tpp[Tpp < Tmin] = Tmin tsnow = Tpp.copy() tsnow[Tpp > Tmax] = Tmax return Tpp, tsnow #BEGIN NASDE MODELS HERE AND BELOW [docs]def susong1999(temperature, precipitation): """ Follows the IPW command mkprecip The precipitation phase, or the amount of precipitation falling as rain or snow, can significantly alter the energy and mass balance of the snowpack, either leading to snow accumulation or inducing melt :cite:`Marks&al:1998` :cite:`Kormos&al:2014`. The precipitation phase and initial snow density are based on the precipitation temperature (the distributed dew point temperature) and are estimated after Susong et al (1999) :cite:`Susong&al:1999`. The table below shows the relationship to precipitation temperature: ========= ======== ============ =============== Min Temp Max Temp Percent snow Snow density [deg C] [deg C] [%] [kg/m^3] ========= ======== ============ =============== -Inf -5 100 75 -5 -3 100 100 -3 -1.5 100 150 -1.5 -0.5 100 175 -0.5 0 75 200 0 0.5 25 250 0.5 Inf 0 0 ========= ======== ============ =============== Args: precipitation - numpy array of precipitation values [mm] temperature - array of temperature values, use dew point temperature if available [degrees C] Returns: dictionary: - **perc_snow** (*numpy.array*) - Percent of the precipitation that is snow in values 0.0-1.0. - **rho_s** (*numpy.array*) - Snow density values in kg/m^3. """ # create a list from the table above t = [] t.append( {'temp_min': -1e309, 'temp_max': -5, 'snow': 1, 'density':75} ) t.append( {'temp_min': -5, 'temp_max': -3, 'snow': 1, 'density':100} ) t.append( {'temp_min': -3, 'temp_max': -1.5, 'snow': 1, 'density':150} ) t.append( {'temp_min': -1.5, 'temp_max': -0.5, 'snow': 1, 'density':175} ) t.append( {'temp_min': -0.5, 'temp_max': 0.0, 'snow': 0.75, 'density':200} ) t.append( {'temp_min': 0.0, 'temp_max': 0.5, 'snow': 0.25, 'density':250} ) t.append( {'temp_min': 0.5, 'temp_max': 1e309, 'snow': 0, 'density':0} ) # preallocate the percent snow (ps) and snow density (sd) ps = np.zeros(precipitation.shape) sd = np.zeros(ps.shape) # if no precipitation return all zeros if np.sum(precipitation) == 0: return ps, sd # determine the indicies and allocate based on the table above for row in t: # get the values between the temperature ranges that have precip ind = ((temperature >= row['temp_min']) & (temperature < row['temp_max'])) # set the percent snow ps[ind] = row['snow'] # set the density sd[ind] = row['density'] # if there is no precipitation at a pixel, don't report a value # this may make isnobal crash, I'm not really sure ps[precipitation == 0] = 0 sd[precipitation == 0] = 0 return {'pcs':ps, 'rho_s':sd} [docs]def piecewise_susong1999(Tpp, precip, Tmax = 0.0, Tmin = -10.0, check_temps=True): """ Follows :func:`~smrf.envphys.snow.susong1999` but is the piecewise form of table shown there. This model adds to the former by accounting for liquid water effect near 0.0 Degrees C. The table was estimated by Danny Marks in 2017 which resulted in the piecewise equations below: Percent Snow: .. math:: \\%_{snow} = \\Bigg \\lbrace{ \\frac{-T}{T_{r0}} P_{cr0} + P_{c0}, \\hfill -0.5^{\\circ} C \\leq T \\leq 0.0^{\\circ} C \\atop \\frac{-T_{pp}}{T_{max} + 1.0} P_{c0} + P_{c0}, \\hfill 0.0^\\circ C \\leq T \\leq T_{max} } Snow Density: .. math:: \\rho_{s} = 50.0 + 1.7 * (T_{pp} + 15.0)^{ex} ex = \\Bigg \\lbrace{ ex_{min} + \\frac{T_{range} + T_{snow} - T_{max}}{T_{range}} * ex_{r}, \\hfill ex < 1.75 \\atop 1.75, \\hfill, ex > 1.75 } Args: Tpp: A numpy array of temperature, use dew point temperature if available [degree C]. precip: A numpy array of precip in millimeters. Tmax: Max temperature that snow density is modeled. Default is 0.0 Degrees C. Tmin: Minimum temperature that snow density is changing. Default is -10.0 Degrees C. check_temps: A boolean value check to apply special temperature constraints, this is done using :func:`~smrf.envphys.snow.check_temperature`. Default is True. Returns: dictionary: - **pcs** (*numpy.array*) - Percent of the precipitation that is snow in values 0.0-1.0. - **rho_s** (*numpy.array*) - Density of the fresh snow in kg/m^3. """ ex_max = 1.75 exr = 0.75 ex_min = 1.0 Tz = 0.0 # again, this shouldn't be needed in this context if check_temps: Tpp, tsnow = check_temperature(Tpp, Tmax=Tmax, Tmin=Tmin) pcs = calc_perc_snow(Tpp,Tmin=Tmin, Tmax = Tmax) # new snow density - no compaction Trange = Tmax - Tmin ex = ex_min + (((Trange + (tsnow - Tmax)) / Trange) * exr) ex[ex > ex_max] = ex_max rho_ns = 50.0 + (1.7 * ((Tpp-Tz) + 15.0)**ex) return {'pcs':pcs, 'rho_s':rho_ns} [docs]def marks2017(Tpp,pp): """ A new accumulated snow density model that accounts for compaction. The model builds upon :func:`~smrf.envphys.snow.piecewise_susong1999` by adding effects from compaction. Of four mechanisms for compaction, this model accounts for compaction by destructive metmorphosism and overburden. These two processes are accounted for by calculating a proportionalility using data from Kojima, Yosida and Mellor. The overburden is in part estimated using total storm deposition, where storms are defined in :func:`~smrf.envphys.storms.tracking_by_station`. Once this is determined the final snow density is applied through the entire storm only varying with hourly temperature. Snow Density: .. math:: \\rho_{s} = \\rho_{ns} + (\\Delta \\rho_{c} + \\Delta \\rho_{m}) \\rho_{ns} Overburden Proportionality: .. math:: \\Delta \\rho_{c} = 0.026 e^{-0.08 (T_{z} - T_{snow})} SWE* e^{-21.0 \\rho_{ns}} Metmorphism Proportionality: .. math:: \\Delta \\rho_{m} = 0.01 c_{11} e^{-0.04 (T_{z} - T_{snow})} c_{11} = c_min + (T_{z} - T_{snow}) C_{factor} + 1.0 Constants: .. math:: C_{factor} = 0.0013 Tz = 0.0 ex_{max} = 1.75 exr = 0.75 ex_{min} = 1.0 c1r = 0.043 c_{min} = 0.0067 c_{fac} = 0.0013 T_{min} = -10.0 T_{max} = 0.0 T_{z} = 0.0 T_{r0} = 0.5 P_{cr0} = 0.25 P_{c0} = 0.75 Args: Tpp: Numpy array of a single hour of temperature, use dew point if available [degrees C]. pp: Numpy array representing the total amount of precip deposited during a storm in millimeters Returns: dictionary: - **rho_s** (*numpy.array*) - Density of the fresh snow in kg/m^3. - **swe** (*numpy.array*) - Snow water equivalent. - **pcs** (*numpy.array*) - Percent of the precipitation that is snow in values 0.0-1.0. - **rho_ns** (*numpy.array*) - Density of the uncompacted snow, which is equivalent to the output from :func:`~smrf.envphys.snow.piecewise_susong1999`. - **d_rho_c** (*numpy.array*) - Prportional coefficient for compaction from overburden. - **d_rho_m** (*numpy.array*) - Proportional coefficient for compaction from melt. - **rho_s** (*numpy.array*) - Final density of the snow [kg/m^3]. - **rho** (*numpy.array*) - Density of the precipitation, which continuously ranges from low density snow to pure liquid water (50-1000 kg/m^3). - **zs** (*numpy.array*) - Snow height added from the precipitation. """ ex_max = 1.75 exr = 0.75 ex_min = 1.0 #c1_min = 0.026 #c1_max = 0.069 c1r = 0.043 c_min = 0.0067 cfac = 0.0013 Tmin = -10.0 Tmax = 0.0 Tz = 0.0 Tr0 = 0.5 Pcr0 = 0.25 Pc0 = 0.75 water = 1000.0 rho_ns = d_rho_m = d_rho_c = zs = rho_s = swe = pcs = np.zeros(Tpp.shape) rho = np.ones(Tpp.shape) if np.any(pp > 0): # check the precipitation temperature Tpp, tsnow = check_temperature(Tpp, Tmax=Tmax, Tmin=Tmin) # Calculate the percent snow and initial uncompacted density result = piecewise_susong1999(Tpp,pp,Tmax=Tmax, Tmin=Tmin) pcs = result['pcs'] rho_orig = result['rho_s'] swe = pp * pcs # Calculate the density only where there is swe swe_ind = swe > 0.0 if np.any(swe_ind): s_pcs = pcs[swe_ind] s_pp = pp[swe_ind] s_swe = swe[swe_ind] s_tpp = Tpp[swe_ind] # transforms to a 1D array, will put back later s_tsnow = tsnow[swe_ind] # transforms to a 1D array, will put back later s_rho_ns = rho_orig[swe_ind] # transforms to a 1D array, will put back later #Convert to a percentage of water s_rho_ns = s_rho_ns/water # proportional total storm mass compaction s_d_rho_c = (0.026 * np.exp(-0.08 * (Tz - s_tsnow)) * s_swe * np.exp(-21.0 * s_rho_ns)) ind = (s_rho_ns * water) >= 100.0 c11 = np.ones(s_rho_ns.shape) #c11[ind] = (c_min + ((Tz - s_tsnow[ind]) * cfac)) + 1.0 c11[ind] = np.exp(-0.046*(s_rho_ns[ind]*1000.0-100.0)) s_d_rho_m = 0.01 * c11 * np.exp(-0.04 * (Tz - s_tsnow)) # compute snow density, depth & combined liquid and snow density s_rho_s = s_rho_ns +((s_d_rho_c + s_d_rho_m) * s_rho_ns) s_zs = s_swe / s_rho_s # a mixture of snow and liquid s_rho = s_rho_s.copy() mix = (s_swe < s_pp) & (s_pcs > 0.0) if np.any(mix): s_rho[mix] = (s_pcs[mix] * s_rho_s[mix]) + (1.0 - s_pcs[mix]) s_rho[s_rho > 1.0] = 1.0 # put the values back into the larger array rho_ns[swe_ind] = s_rho_ns d_rho_m[swe_ind] = s_d_rho_m d_rho_c[swe_ind] = s_d_rho_c zs[swe_ind] = s_zs rho_s[swe_ind]= s_rho_s rho[swe_ind] = s_rho pcs[swe_ind] = s_pcs # convert densities from proportions, to kg/m^3 or mm/m^2 rho_ns = rho_ns * water rho_s = rho_s * water rho = rho * water return {'swe':swe, 'pcs':pcs,'rho_ns': rho_ns, 'd_rho_c' : d_rho_c, 'd_rho_m' : d_rho_m, 'rho_s' : rho_s, 'rho':rho, 'zs':zs} #Available NASDE Models should be put in this dictionary as well: available_models = {"susong1999":susong1999, "piecewise_susong1999": piecewise_susong1999, "marks2017":marks2017} if __name__ == '__main__': print("\nNothing implemented here.")