Merge pull request #776 from dgergel/feature/fix_deltacc_bug

Feature/fix deltacc bug

docs/Development/ReleaseNotes.md

@@ -115,6 +115,10 @@ To check which release of VIC you are running:
 
     Previously this was set in the surface_fluxes.c numerics routine for ground-canopy iterations, which meant that that routine had to be altered to change the maximum number of iterations. It has now been moved to the parameters struct so that it can be overriden in the constants file. 
 
+10. Updated new snow density function by adding a cap to new snow density that is set in the parameters struct by the parameter SNOW_NEW_SNOW_DENS_MAX ([GH#776](https://github.com/UW-Hydro/VIC/pull/776))
+
+    Previously the change in cold content of the snowpack term (deltaCC in the snow_data_struct) would get unreasonably large if the Hedstrom and Pomeroy 1998 equation used to calculate snow density, which depends only on air temperature, was calculated with air temperatures above about 2 deg C. We use this term to calculate the ground flux from the snowpack and snow depth, which resulted in extremely small snow depths and unreasonably large ground fluxes from the snowpack (and thus changes in snowpack cold content). Now there is a cap on new snow density with the new parameter SNOW_NEW_SNOW_DENS_MAX as well as a snow depth below which we disregard the ground flux from the snowpack (1.e-8). 
+
 10. Miscellaneous clean-up:
 
     [GH#723](https://github.com/UW-Hydro/VIC/pull/723)

docs/Documentation/Constants.md

@@ -111,6 +111,8 @@ The table below lists the constants available for manipulation via the `CONSTANT
 | SNOW_MAX_SURFACE_SWE         |             |
 | SNOW_LIQUID_WATER_CAPACITY   |             |
 | SNOW_NEW_SNOW_DENSITY        |             |
+| SNOW_NEW_SNOW_DENS_MAX       |             |
+| SNOW_DEPTH_THRES             |             |
 | SNOW_DENS_DMLIMIT            |             |
 | SNOW_DENS_DMLIMIT_FACTOR     |             |
 | SNOW_DENS_MAX_CHANGE         |             |

vic/drivers/cesm/bld/vic.constants.txt

@@ -122,6 +122,8 @@
 # SNOW_MAX_SURFACE_SWE 0.125
 # SNOW_LIQUID_WATER_CAPACITY 0.035
 # SNOW_NEW_SNOW_DENSITY 50.0
+# SNOW_NEW_SNOW_DENS_MAX 400.0
+# SNOW_DEPTH_THRES 1.e-8
 # SNOW_DENS_DMLIMIT 100.0
 # SNOW_DENS_MAX_CHANGE 0.9
 # SNOW_DENS_ETA0 3.6e6

vic/drivers/shared_all/src/get_parameters.c

@@ -326,6 +326,12 @@ get_parameters(FILE *paramfile)
             else if (strcasecmp("SNOW_NEW_SNOW_DENSITY", optstr) == 0) {
                 sscanf(cmdstr, "%*s %lf", &param.SNOW_NEW_SNOW_DENSITY);
             }
+            else if (strcasecmp("SNOW_NEW_SNOW_DENS_MAX", optstr) == 0) {
+                sscanf(cmdstr, "%*s %lf", &param.SNOW_NEW_SNOW_DENS_MAX);
+            }
+            else if (strcasecmp("SNOW_NEW_SNOW_DENS_MAX", optstr) == 0) {
+                sscanf(cmdstr, "%*s %lf", &param.SNOW_NEW_SNOW_DENS_MAX);
+            }
             else if (strcasecmp("SNOW_DENS_DMLIMIT", optstr) == 0) {
                 sscanf(cmdstr, "%*s %lf", &param.SNOW_DENS_DMLIMIT);
             }
@@ -758,10 +764,19 @@ validate_parameters()
         log_err(
             "SNOW_NEW_SNOW_DENSITY must be defined on the interval [0, inf) (kg/m^3)");
     }
+    if (!(param.SNOW_DEPTH_THRES >= 0.)) {
+        log_err(
+            "SNOW_DEPTH_THRES must be defined on the interval [0, inf) (m)");
+    }
     if (!(param.SNOW_DENS_DMLIMIT >= 0.)) {
         log_err(
             "SNOW_DENS_DMLIMIT must be defined on the interval [0, inf) (kg/m^3)");
     }
+    if (!(param.SNOW_NEW_SNOW_DENS_MAX >= 0. &&
+          param.SNOW_NEW_SNOW_DENS_MAX <= 700.)) {
+        log_err(
+            "SNOW_NEW_SNOW_DENS_MAX must be defined on the interval [0, 700) (kg/m^3)");
+    }
     if (!(param.SNOW_DENS_MAX_CHANGE >= 0 && param.SNOW_DENS_MAX_CHANGE <= 1)) {
         log_err("SNOW_DENS_MAX_CHANGE must be defined on the interval [0,1] (-)")
     }

vic/drivers/shared_all/src/initialize_parameters.c

@@ -151,6 +151,8 @@ initialize_parameters()
     param.SNOW_MAX_SURFACE_SWE = 0.125;
     param.SNOW_LIQUID_WATER_CAPACITY = 0.035;
     param.SNOW_NEW_SNOW_DENSITY = 50.0;
+    param.SNOW_NEW_SNOW_DENS_MAX = 400.0;
+    param.SNOW_DEPTH_THRES = 1.e-8;
     param.SNOW_DENS_DMLIMIT = 100.0;
     param.SNOW_DENS_DMLIMIT_FACTOR = 1.15;
     param.SNOW_DENS_MAX_CHANGE = 0.9;

vic/drivers/shared_all/src/print_library_shared.c

@@ -753,6 +753,10 @@ print_parameters(parameters_struct *param)
             param->SNOW_LIQUID_WATER_CAPACITY);
     fprintf(LOG_DEST, "\tSNOW_NEW_SNOW_DENSITY: %.4f\n",
             param->SNOW_NEW_SNOW_DENSITY);
+    fprintf(LOG_DEST, "\tSNOW_NEW_SNOW_DENS_MAX: %.4f\n",
+            param->SNOW_NEW_SNOW_DENS_MAX);
+    fprintf(LOG_DEST, "\tSNOW_DEPTH_THRES: %.12f\n",
+            param->SNOW_DEPTH_THRES);
     fprintf(LOG_DEST, "\tSNOW_DENS_DMLIMIT: %.4f\n", param->SNOW_DENS_DMLIMIT);
     fprintf(LOG_DEST, "\tSNOW_DENS_MAX_CHANGE: %.4f\n",
             param->SNOW_DENS_MAX_CHANGE);

vic/drivers/shared_image/src/vic_mpi_support.c

@@ -755,7 +755,7 @@ create_MPI_param_struct_type(MPI_Datatype *mpi_type)
     MPI_Datatype   *mpi_types;
 
     // nitems has to equal the number of elements in parameters_struct
-    nitems = 154;
+    nitems = 156;
     blocklengths = malloc(nitems * sizeof(*blocklengths));
     check_alloc_status(blocklengths, "Memory allocation error.");
 
@@ -1133,6 +1133,14 @@ create_MPI_param_struct_type(MPI_Datatype *mpi_type)
     offsets[i] = offsetof(parameters_struct, SNOW_NEW_SNOW_DENSITY);
     mpi_types[i++] = MPI_DOUBLE;
 
+    // double SNOW_NEW_SNOW_DENS_MAX
+    offsets[i] = offsetof(parameters_struct, SNOW_NEW_SNOW_DENS_MAX);
+    mpi_types[i++] = MPI_DOUBLE;
+
+    // double SNOW_DEPTH_THRES
+    offsets[i] = offsetof(parameters_struct, SNOW_DEPTH_THRES);
+    mpi_types[i++] = MPI_DOUBLE;
+
     // double SNOW_DENS_DMLIMIT
     offsets[i] = offsetof(parameters_struct, SNOW_DENS_DMLIMIT);
     mpi_types[i++] = MPI_DOUBLE;

vic/vic_run/include/vic_def.h

@@ -453,6 +453,8 @@ typedef struct {
     double SNOW_MAX_SURFACE_SWE;  /**< maximum depth of the surface layer in water equivalent (m) */
     double SNOW_LIQUID_WATER_CAPACITY;  /**< water holding capacity of snow as a fraction of snow-water-equivalent */
     double SNOW_NEW_SNOW_DENSITY;  /**< density of new fallen snow */
+    double SNOW_NEW_SNOW_DENS_MAX; /**< new snow density max for Hedstrom and Pomeroy 1998 equation [Warren et al. 1999, Bormann et al. 2013, Maidment Figure 7.2.3] */
+    double SNOW_DEPTH_THRES;  /**< Snow depth threshold below which we do not consider the ground flux out of the snowpack in calculating change in cold content (m) */
     double SNOW_DENS_DMLIMIT;  /**< Density limit used in calculation of destructive metamorphism (kg/m^3) */
     double SNOW_DENS_DMLIMIT_FACTOR;  /**< Density limit factor used in calculation of destructive metamorphism (kg/m^3) */
     double SNOW_DENS_MAX_CHANGE;  /**< maximum change in snowfall depth (fraction of swe) */

vic/vic_run/src/SnowPackEnergyBalance.c

@@ -3,6 +3,10 @@
  *
  * Calculate snow pack energy balance
  *
+ * Based on the SnowPackEnergyBalance function in DHSVM
+ * Reference: Bras, R.A., Hydrology, an introduction to hydrologic science,
+ * Addison Wesley, Inc., Reading, etc., 1990.
+ *
  * @section LICENSE
  *
  * The Variable Infiltration Capacity (VIC) macroscale hydrological model
@@ -238,20 +242,20 @@ SnowPackEnergyBalance(double  TSurf,
        Equation 7.3.12 from H.B.H. for rain falling on melting snowpack */
 
     if (TMean == 0.) {
-        *AdvectedEnergy = (CONST_CPFW * CONST_RHOFW * (Tair) * Rain) / (Dt);
+        *AdvectedEnergy = (CONST_CPFW * CONST_RHOFW * (Tair) * Rain) / Dt;
     }
     else {
         *AdvectedEnergy = 0.;
     }
 
     /* Calculate change in cold content */
     *DeltaColdContent = CONST_VCPICE_WQ * SweSurfaceLayer *
-                        (TSurf - OldTSurf) / (Dt);
+                        (TSurf - OldTSurf) / Dt;
 
     /* Calculate Ground Heat Flux */
-    if (SnowDepth > 0.) {
-        *GroundFlux = 2.9302e-6 * SnowDensity * SnowDensity *
-                      (TGrnd - TMean) / SnowDepth / (Dt);
+    if (SnowDepth > param.SNOW_DEPTH_THRES) {
+        *GroundFlux = param.SNOW_CONDUCT * pow(SnowDensity, 2.) *
+                      (TGrnd - TMean) / SnowDepth / Dt;
     }
     else {
         *GroundFlux = 0;
@@ -263,7 +267,7 @@ SnowPackEnergyBalance(double  TSurf,
                *AdvectedEnergy + *GroundFlux - *DeltaColdContent +
                *AdvectedSensibleHeat;
 
-    *RefreezeEnergy = (SurfaceLiquidWater * CONST_LATICE * Density) / (Dt);
+    *RefreezeEnergy = (SurfaceLiquidWater * CONST_LATICE * Density) / Dt;
 
     if (TSurf == 0.0 && RestTerm > -(*RefreezeEnergy)) {
         *RefreezeEnergy = -RestTerm; /* available energy input over cold content

vic/vic_run/src/snow_utility.c

@@ -235,6 +235,12 @@ new_snow_density(double air_temp)
         log_err("Unknown SNOW_DENSITY option");
     }
 
+    // cap new snow density to prevent the calculation from
+    // becoming unphysical
+    if (density_new > param.SNOW_NEW_SNOW_DENS_MAX) {
+        density_new = param.SNOW_NEW_SNOW_DENS_MAX;
+    }
+
     return (density_new);
 }