snowpack-model / snowpack Goto Github PK

I'm noticing some strange differences between WINDEROSIONDEPOSITION and delta SWE.

For example, after a few time steps, the time sum of WINDEROSIONDEPOSITION looks like this. Sensible right?

But then delta SWE looks like this.

Somehow I think the different variables in SnowpackInterface.cc are not being updated properly.

Top row of netcdf is not written with data, but contains nodata. This can either be that the field is shifted one column down, or that the top row is never filled with data.

As discussed in Amory et al. (2017), drag coefficients, and thus roughness length can be considered a function of air temperature, to separate between summer and winter. See particularly eq. 7 for the temperature dependence, and eq. 4 to convert drag coefficients in roughness length, under the assumption of a neutral boundary.

Citation: Amory, Charles; Gallée, Hubert; Naaim-Bouvet, Florence; Favier, Vincent; Vignon, Etienne; Picard, Ghislain; Trouvilliez, Alexandre; Piard, Luc; Genthon, Christophe; Bellot, Hervé (2017). Seasonal Variations in Drag Coefficient over a Sastrugi-Covered Snowfield in Coastal East Antarctica. Boundary-Layer Meteorology, 164(1), 107–133. doi:10.1007/s10546-017-0242-5

Delta SWE is identical to (MS_HNW - MS_WIND - ET - MS_SNOWPACK_RUNOFF). But is not equal to (pr - WINDEROSIONDEPOSITION - ET - MS_SNOWPACK_RUNOFF)

Based off my understanding of the variables I would expect both expressions to be equal to detla SWE. Do you agree? If so, could the discrepancy be due to a programing error? My impression is that numerical errors arising from explicit snow drift should be encapsulated in delta swe and WINDEROSIONDEPOSITION, and therefore should not explain the difference.

The dimensionality of two parameters (0.0014 and 205 below) are currently incorrect. Magnitudes currently reflect units of g/cm/s, but should be transformed to units of kg/m/s.

Additional context from reviewer comment:

l.124 (eq.4): The numerator of this equation corresponds to the expression proposed by Sørensen (1991) for the transport rate (see page 75 of the article, eq. 3.22). This expression is in units of g/cm/s (this is stated at the end of page 72 of the article, below equation 3.9). This poses a problem because the coefficients 0.0014 and 205 are not dimensionless values - 205 has velocity units (cm/s) and 0.0014 has units of s2/cm. If we want to express Q in units of kg/m/s, these factors should change to 2.05 and 0.14, respectively. The dimensionally correct expression predicts much higher values of Q and its validity to model snow saltation is still to be assessed. This issue with the Sørensen's expression was previously pointed out in the PhD thesis of Vionnet (2012), page 103, Fig. 5.3 (french only). In the mentioned PhD thesis as well as in Vionnet et al. (2014), the use of the latest expression of Sørensen (2004) is proposed. This can be a good option for Alpine3D as it does not deviate significantly from the dimensionally wrong Sørensen equation (see Fig.5.3 in the PhD thesis). Independently of the approach chosen by the authors, I believe it is advisable to present Q - the numerator of eq.4 - in a separate equation and cite the respective article.

Discussion of potential solutions in Vincent Vionnet's 2012 PhD thesis (page 103, Fig 5.3)

Alpine3D is occasionally considering a pixel a glacier, and then forcing albedo to 0.3. The snow layers don't indicate ice/glacier (i.e., mk 7, 8 or snowtype 880).

This issue doesn't seem to be the case for the DML simulation.

We should add the following time-invariant fields to Alpine-3D netCDF output:

1. Latitude/longitude of each grid cell in WGS-84 (in addition to X and Y coordinates in EPSG 3031 which Alpine-3D already writes)
2. Grid cell area (m^2)
3. Grid cell elevation (m)

netCDF metadata indicate SWE units are meters liquid water equivalent, but values are around 1000. So I guess they are in mm?

xarray.DataArray'swe'northing: 107easting: 168
1040.2449 1040.1312 1039.4742 1039.2972 1039.1523 ... nan nan nan nan
Coordinates:
time
()
datetime64[ns]
1980-01-09T12:00:00
northing
(northing)
float32
-109000.0 -108000.0 ... -3000.0
easting
(easting)
float32
-1542000.0 ... -1375000.0
Attributes:
standard_name :
lwe_thickness_of_surface_snow_amount
units :
m

We should add the time variable fields listed here to Alpine-3D netCDF output. Examples include:

1. Top 1 m snow density
2. Liquid water content
3. Snow height change

In SnowpackInterface.cc lines 1538-1539 (divergenceSnowDrift branch):

Q_x(ix, iy) = (tmp_ErodedMass(ix, iy) * L * U(ix, iy)) / (dt * grid_VW(ix, iy));
Q_y(ix, iy) = (tmp_ErodedMass(ix, iy) * L * V(ix, iy)) / (dt * grid_VW(ix, iy));

Q_x and Q_y are both very small, ~5e-4. But C++ stores the values as zero.

How can I increase the numerical precision? I assume Q_x and Q_y should be doubles? Not sure what they are initialized as

mio::Grid2DObject Q_x( ErodedMass, 0. ); // Mass flux through a vertical gate in the x direction (kg/m/s)
mio::Grid2DObject Q_y( ErodedMass, 0. ); // Mass flux through a vertical gate in the y direction (kg/m/s) 

I'm working on adding surface mass balance (SMB) to SnowpackInterfaceWorker.cc

To do so I am adding a new case to void SnowpackInterfaceWorker::fillGrids(const size_t& ii, const size_t& jj, const CurrentMeteo& meteoPixel, const SnowStation& snowPixel, const SurfaceFluxes& surfaceFlux). Does this look roughly right? I'm not farmiliar with tests for negative numbers, float/int, etc in C++. So any help/review would be appreciated!

So far I've tried to account for precipitation, sublimation, and meltwater runoff. I don't quite understand how to add drifting snow erosion... Am I correct that at the moment, SNOWPACK has no drifting snow sublimation term?

case SnGrids::SMB:
   value = meteoPixel.psum  - (surfaceFlux.mass[SurfaceFluxes::MS_SUBLIMATION] / snowPixel.cos_sl) - surfaceFlux.mass[SurfaceFluxes::MS_SNOWPACK_RUNOFF]; break;

Add a field for the age of a layer to the *pro file.

Add feedback between HS and DEM as an option to Alpine3D.

The optimized SNOWPACK / Alpine3D (which aims to reuse classes) doesn't work properly. Some weird simulation results show up. No crashes, and valgrind doesn't complain.

First, it should be tried to add copy-constructor and assignment operator.

Quick fix for the described issue

diff --git a/Source/meteoio/meteoio/dataClasses/DEMObject.cc b/Source/meteoio/meteoio/dataClasses/DEMObject.cc
index 2c3cf2a..f8eb452 100644
--- a/Source/meteoio/meteoio/dataClasses/DEMObject.cc
+++ b/Source/meteoio/meteoio/dataClasses/DEMObject.cc
@@ -431,6 +431,8 @@ void DEMObject::CalculateAziSlopeCurve(slope_type algorithm) {
 //This computes the slope and the aspect at a given cell as well as the x and y components of the normal vector
 	double A[4][4]; //table to store neigbouring heights: 3x3 matrix but we want to start at [1][1]
 	                //we use matrix notation: A[y][x]
+	double A2[4][4]; //table to store neigbouring heights: 3x3 matrix but we want to start at [1][1]
+	                //we use matrix notation: A[y][x]
 	if (algorithm==DFLT) {
 		algorithm = dflt_algorithm;
 	}
@@ -467,10 +469,11 @@ void DEMObject::CalculateAziSlopeCurve(slope_type algorithm) {
 				}
 			} else {
 				getNeighbours(i, j, A);
+				getNeighbours2(i, j, A2);
 				double new_slope, new_Nx, new_Ny, new_Nz;
 				(this->*CalculateSlope)(A, new_slope, new_Nx, new_Ny, new_Nz);
 				const double new_azi = CalculateAzimuth(new_Nx, new_Ny, new_Nz, new_slope);
-				const double new_curvature = getCurvature(A);
+				const double new_curvature = getCurvature(A2);
 				if (update_flag&SLOPE) {
 					slope(i,j) = new_slope;
 					azi(i,j) = new_azi;
@@ -782,6 +785,34 @@ void DEMObject::getNeighbours(const size_t& i, const size_t& j, double A[4][4])
 	}
 }
 
+void DEMObject::getNeighbours2(const size_t& i, const size_t& j, double A[4][4]) const
+{ //this fills a 3x3 table containing the neighboring values
+	const int n_radius = 10; // Number of grid cells which contain the topographic slope curvature length scale.
+
+	if ((i>n_radius-1 && i<(getNx()-n_radius)) && (j>n_radius-1 && j<(getNy()-n_radius))) {
+		//this is the normal case
+		A[1][1] = grid2D(i-n_radius, j+n_radius);
+		A[1][2] = grid2D(i, j+n_radius);
+		A[1][3] = grid2D(i+n_radius, j+n_radius);
+		A[2][1] = grid2D(i-n_radius, j);
+		A[2][2] = grid2D(i, j);
+		A[2][3] = grid2D(i+n_radius, j);
+		A[3][1] = grid2D(i-n_radius, j-n_radius);
+		A[3][2] = grid2D(i, j-n_radius);
+		A[3][3] = grid2D(i+n_radius, j-n_radius);
+	} else { //somewhere on the border
+		A[1][1] = safeGet((signed)i-n_radius, (signed)j+n_radius);
+		A[1][2] = safeGet((signed)i, (signed)j+n_radius);
+		A[1][3] = safeGet((signed)i+n_radius, (signed)j+n_radius);
+		A[2][1] = safeGet((signed)i-n_radius, (signed)j);
+		A[2][2] = safeGet((signed)i, (signed)j);
+		A[2][3] = safeGet((signed)i+n_radius, (signed)j);
+		A[3][1] = safeGet((signed)i-n_radius, (signed)j-n_radius);
+		A[3][2] = safeGet((signed)i, (signed)j-n_radius);
+		A[3][3] = safeGet((signed)i+n_radius, (signed)j-n_radius);
+	}
+}
+
 double DEMObject::safeGet(const int& i, const int& j) const
 {//this function would allow reading the value of *any* point,
 //that is, even for coordinates outside of the grid (where it would return nodata)
diff --git a/Source/meteoio/meteoio/dataClasses/DEMObject.h b/Source/meteoio/meteoio/dataClasses/DEMObject.h
index 3d1fd31..0e1c298 100644
--- a/Source/meteoio/meteoio/dataClasses/DEMObject.h
+++ b/Source/meteoio/meteoio/dataClasses/DEMObject.h
@@ -146,6 +146,7 @@ class DEMObject : public Grid2DObject {
 		static void surfaceGradient(double& dx_sum, double& dy_sum, double A[4][4]);
 		static double avgHeight(const double& z1, const double &z2, const double& z3);
 		void getNeighbours(const size_t& i, const size_t& j, double A[4][4]) const;
+		void getNeighbours2(const size_t& i, const size_t& j, double A[4][4]) const;
 		double safeGet(const int& i, const int& j) const;
 
 		double max_shade_distance; ///< maximum distance to look for when computing a horizon

We should write a consistent wind direction function, instead of writing it in several places. Here is a start

// Calculate wind direction from u and v components of wind.
// Wind direction is calculated in degrees in the meteorological sense. For example, a wind coming from the north has a wind direction of zero degrees. 
// u and v components denote the horizontal and vertical components respectively and indicate the direction the wind is going. I.e. a wind coming from the north would have u = 0 and v < 0.
wind_direction = fmod( atan2(u, v) * Cst::to_deg + 180., 360.);

In the driftingsnow branch, we establish an adaptive time step for the explicit snow drift scheme in order to simultaneously minimize computational expense and preserve numerical stability.

Currently we have implemented a particle slow down term. However, this has not yet been accounted for when calculating the sub time step. In the case of a constant slowdown, this could be fixed with

double sub_dt = std::min(C_max * dx / (sqrt(2.) * grid_VW.grid2D.getMax() - particle_slowdown), dt); // Sub time step

or in the case of a multiplicative slowdown

double sub_dt = std::min(C_max * dx / (particle_slowdown * sqrt(2.) * grid_VW.grid2D.getMax()), dt); // Sub time step

Add support for MS_WIND in Alpine3D.

Currently, all specified grids are written to a single NetCDF file. For my purposes, these files are becoming too large. For me, it would be preferred to write individual files for each individual year. We could also consider adding an experiment id with information about forcing, spinup, etc.

e.g.

WAIS_MERRA2_100m_spinup_T2M_1980_hourly.nc 
WAIS_MERRA2_100m_spinup_SWE_1980_hourly.nc
....
....
WAIS_MERRA2_100m_spinup_T2M_2020_hourly.nc 
WAIS_MERRA2_100m_spinup_SWE_2020_hourly.nc

To make analysis easier, it would be nice to have FAC and SMB built into SNOWPACK output. This way we don't have to calculate it ourselves (and possibly screw it up)!

SMB would look something like this (right?)

smet['SMB'] = smet['MS_Snow'] - smet['MS_Wind'] + smet['rainfall'] + smet['sublimation'] + smet['evap'] + smet['runoff']

Firn air content could be calculated by vertically integrating element air content.