Skip to content

Commit

Permalink
update snow
Browse files Browse the repository at this point in the history
  • Loading branch information
@ftomei
ftomei committed Dec 12, 2023
1 parent 44a2b1f commit 32182c7
Show file tree
Hide file tree
Showing 8 changed files with 130 additions and 48 deletions.
1 change: 1 addition & 0 deletions DATA/PROJECT/Troy/SETTINGS/parameters.ini
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,3 +73,4 @@ snowSkinThickness=0.003
snowVegetationHeight=0.1
soilAlbedo=0.2
skinThickness=0.02
snowSurfaceDampingDepth=0.05
104 changes: 74 additions & 30 deletions agrolib/snow/snow.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,13 +40,13 @@ Crit3DSnowParameters::Crit3DSnowParameters()
void Crit3DSnowParameters::initialize()
{
// default values
skinThickness = 0.02; /*!< [m] */ // LC: VERSIONI DIVERSE IN BROOKS: 3mm (nel testo), 2-3cm (nel codice)
skinThickness = 0.02; /*!< [m] */ // VERSIONI DIVERSE IN BROOKS: 3mm (nel testo), 2-3cm (nel codice)
soilAlbedo = 0.2; /*!< [-] bare soil */
snowVegetationHeight = 1;
snowWaterHoldingCapacity = 0.05;
snowMaxWaterContent = 0.1;
tempMaxWithSnow = 2;
tempMinWithRain = -0.5;
snowVegetationHeight = 1; /*!< [m] */
snowWaterHoldingCapacity = 0.05; /*!< [-] */
tempMaxWithSnow = 2; /*!< [°C] */
tempMinWithRain = -0.5; /*!< [°C] */
snowSurfaceDampingDepth = 0.05; /*!< [m] */ // VERSIONI DIVERSE IN BROOKS: 0.05 - 0.15
}


Expand Down Expand Up @@ -145,7 +145,6 @@ void Crit3DSnow::computeSnowBrooksModel()
{
double solarRadTot;
double cloudCover; /*!< [-] */
double prevIceContent, prevLWaterContent; /*!< [mm] */
double currentRatio;

double airActualVapDensity; /*!< [kg m-3] */
Expand All @@ -171,7 +170,7 @@ void Crit3DSnow::computeSnowBrooksModel()

// free water
bool isWater = false;
if ((_surfaceWaterContent / 1000) > snowParameters.snowMaxWaterContent ) /*!< [m] */
if (_surfaceWaterContent > 100. ) /*!< [mm] acqua libera (fiumi - torrenti) */
isWater = true;

if (isWater || (! checkValidPoint()))
Expand Down Expand Up @@ -216,32 +215,31 @@ void Crit3DSnow::computeSnowBrooksModel()
double prevInternalEnergy = _internalEnergy;
double prevSurfaceEnergy = _surfaceEnergy;
double prevSurfaceTemp = _surfaceTemp;
double prevIceContent = _iceContent;
double prevLWaterContent = _liquidWaterContent;

/*--------------------------------------------------------------------
// controlli di coerenza per eventuali modifiche manuali su mappa SWE
// -------------------------------------------------------------------*/
if (previousSWE > 0)
{
prevIceContent = _iceContent;
prevLWaterContent = _liquidWaterContent;

if ( (prevIceContent <= 0) && (prevLWaterContent <= 0) )
/*--------------------------------------------------------------------
// controlli di coerenza per eventuali modifiche manuali su mappa SWE
// -------------------------------------------------------------------*/
if (prevIceContent <= 0 && prevLWaterContent <= 0)
{
// neve aggiunta
prevIceContent = previousSWE;
prevLWaterContent = previousSWE * snowParameters.snowWaterHoldingCapacity / (1 - snowParameters.snowWaterHoldingCapacity);

// Pag. 53 formula 3.23
prevInternalEnergy = -(previousSWE / 1000) * LATENT_HEAT_FUSION * WATER_DENSITY;
prevInternalEnergy = -previousSWE * 0.001 * LATENT_HEAT_FUSION * WATER_DENSITY;

// stato fittizio:: neve recente prossima alla fusione, con una settimana di età
_ageOfSnow = 7;
prevSurfaceTemp = std::min(prevSurfaceTemp, 0.);
prevSurfaceEnergy = std::min(prevSurfaceEnergy, 0.);
prevSurfaceEnergy = computeSurfaceEnergySnow(prevSurfaceTemp, std::min(previousSWE, snowParameters.skinThickness));

_ageOfSnow = 1;
}

/*! check on sum */
currentRatio = previousSWE / (prevIceContent + prevLWaterContent);
if (fabs(currentRatio - 1) > 0.001)
if (! isEqual(currentRatio, 1) )
{
prevIceContent = prevIceContent * currentRatio;
prevLWaterContent = prevLWaterContent * currentRatio;
Expand All @@ -254,7 +252,31 @@ void Crit3DSnow::computeSnowBrooksModel()
_ageOfSnow = NODATA;
}

/*! \brief Vapor Density and Roughness Calculations */
/*! \brief check on soil internal energy - added by ftomei */

if ( previousSWE < EPSILON )
{
// no snow condition: all soil
double estInternalEnergy = computeInternalEnergySoil(prevSurfaceTemp, DEFAULT_BULK_DENSITY); // [kJ m-2]
double absDifference = fabs(estInternalEnergy - prevInternalEnergy); // [kJ m-2]

// the difference is very high
if (absDifference > 1000)
{
// check to avoid division by zero
if ( isEqual(estInternalEnergy, 0) )
{
estInternalEnergy = EPSILON;
}
double ratio = prevInternalEnergy / estInternalEnergy; // [-]
if ( (ratio < 0.5) || (ratio > 2))
{
prevInternalEnergy = (prevInternalEnergy + estInternalEnergy) * 0.5;
}
}
}

/*! \brief compute Roughness and Vapor Density */

// brooks originale
if ( previousSWE > SNOW_MINIMUM_HEIGHT)
Expand Down Expand Up @@ -449,16 +471,25 @@ void Crit3DSnow::computeSnowBrooksModel()
/*! Snow water equivalent */
_snowWaterEquivalent = _iceContent + _liquidWaterContent;

/*! surface energy [kJ m-2] and surface temeprature [°C] */
// all snow
double surfaceEnergySnow = std::min(0., prevSurfaceEnergy + (QTotal + Qr) * (snowParameters.skinThickness / SNOW_DAMPING_DEPTH) );
/*! surface energy [kJ m-2] and surface temperature [°C] */
double surfaceEnergySnow;
// snow
if (_snowWaterEquivalent > 0 && fabs(_internalEnergy) < EPSILON)
{
surfaceEnergySnow = 0.;
}
else
{
double snowRatio = std::min(snowWaterEquivalent * 0.001, snowParameters.skinThickness) / snowParameters.snowSurfaceDampingDepth;
surfaceEnergySnow = std::min(0., prevSurfaceEnergy + (QTotal + Qr) * snowRatio);
}
double surfaceTempSnow = surfaceEnergySnow / (WATER_DENSITY * SNOW_SPECIFIC_HEAT * snowParameters.skinThickness);

// all soil
double surfaceEnergySoil = prevSurfaceEnergy + (QTotal + Qr) * (snowParameters.skinThickness / SOIL_DAMPING_DEPTH);
double surfaceTempSoil = surfaceEnergySoil / (DEFAULT_BULK_DENSITY * SOIL_SPECIFIC_HEAT * snowParameters.skinThickness);

double snowDepthRatio = 5.;
double snowDepthRatio = 4.;
double snowFraction = std::min(_snowWaterEquivalent * snowDepthRatio / 1000., snowParameters.skinThickness) / snowParameters.skinThickness;

_surfaceEnergy = (surfaceEnergySnow * snowFraction) + surfaceEnergySoil * (1 - snowFraction);
Expand Down Expand Up @@ -611,14 +642,27 @@ double aerodynamicResistanceCampbell77(bool isSnow , double zRefWind, double win


// pag. 54 (3.29)
double computeInternalEnergy(double initSoilPackTemp,int bulkDensity, double initSWE)
double computeInternalEnergy(double soilTemperature, int bulkDensity, double swe)
{
return soilTemperature * ( WATER_DENSITY * SNOW_SPECIFIC_HEAT * swe * 0.001
+ bulkDensity * SOIL_SPECIFIC_HEAT * SOIL_DAMPING_DEPTH );
}


double computeInternalEnergySoil(double soilTemperature, int bulkDensity)
{
return soilTemperature * bulkDensity * SOIL_SPECIFIC_HEAT * SOIL_DAMPING_DEPTH;
}


double computeSurfaceEnergySnow(double surfaceTemperature, double skinThickness)
{
return initSoilPackTemp * (WATER_DENSITY * SNOW_SPECIFIC_HEAT * initSWE + bulkDensity * SOIL_SPECIFIC_HEAT * SNOW_DAMPING_DEPTH);
return surfaceTemperature * WATER_DENSITY * SNOW_SPECIFIC_HEAT * skinThickness;
}


double computeSurfaceEnergy(double initSnowSurfaceTemp, double skinThickness)
double computeSurfaceEnergySoil(double surfaceTemperature, double skinThickness)
{
return initSnowSurfaceTemp * WATER_DENSITY * SNOW_SPECIFIC_HEAT * skinThickness;
return surfaceTemperature * DEFAULT_BULK_DENSITY * SOIL_SPECIFIC_HEAT * skinThickness;
}

11 changes: 7 additions & 4 deletions agrolib/snow/snow.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@
#define SOIL_SPECIFIC_HEAT 1.4 /*!< [KJ kg-1 °C-1] wet soil */
#define DEFAULT_BULK_DENSITY 1350 /*!< [kg m-3] */
#define SOIL_DAMPING_DEPTH 0.3 /*!< [m] */
#define SNOW_DAMPING_DEPTH 0.05 /*!< [m] */
#define SNOW_MINIMUM_HEIGHT 1. /*!< [mm] */


Expand All @@ -33,9 +32,9 @@
double soilAlbedo; /*!< [-] bare soil */
double snowVegetationHeight; /*!< [m] height of vegetation */
double snowWaterHoldingCapacity; /*!< [-] percentuale di acqua libera che il manto nevoso può trattenere */
double snowMaxWaterContent; /*!< [m] acqua libera (torrenti, laghetti) */
double tempMaxWithSnow; /*!< [°C] */
double tempMinWithRain; /*!< [°C] */
double snowSurfaceDampingDepth; /*!< [m] */

Crit3DSnowParameters();

Expand Down Expand Up @@ -112,8 +111,12 @@


double aerodynamicResistanceCampbell77(bool isSnow , double zRefWind, double windSpeed, double vegetativeHeight);
double computeInternalEnergy(double initSoilPackTemp,int bulkDensity, double initSWE);
double computeSurfaceEnergy(double initSnowSurfaceTemp, double skinThickness);

double computeInternalEnergy(double soilTemperature,int bulkDensity, double swe);
double computeInternalEnergySoil(double soilTemperature, int bulkDensity);

double computeSurfaceEnergySnow(double surfaceTemperature, double skinThickness);
double computeSurfaceEnergySoil(double surfaceTemperature, double skinThickness);


#endif // SNOW_H
Expand Down
9 changes: 8 additions & 1 deletion agrolib/snow/snowMaps.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -172,7 +172,14 @@ void Crit3DSnowMaps::resetSnowModel(double skinThickness)

_snowSurfaceTempMap->value[row][col] = float(_initSnowSurfaceTemp);

_surfaceEnergyMap->value[row][col] = float(computeSurfaceEnergy(_initSnowSurfaceTemp, skinThickness));
if (initSWE > 0)
{
_surfaceEnergyMap->value[row][col] = float(computeSurfaceEnergySnow(_initSnowSurfaceTemp, skinThickness));
}
else
{
_surfaceEnergyMap->value[row][col] = float(computeSurfaceEnergySoil(_initSnowSurfaceTemp, skinThickness));
}

_internalEnergyMap->value[row][col] = float(computeInternalEnergy(_initSoilPackTemp, surfaceBulkDensity, initSWE/1000.));

Expand Down
6 changes: 6 additions & 0 deletions bin/CRITERIA3D/criteria3DProject.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -512,12 +512,17 @@ bool Crit3DProject::loadCriteria3DParameters()
{
snowModel.snowParameters.soilAlbedo = parameters->value("soilAlbedo").toDouble();
}
if (parameters->contains("snowSurfaceDampingDepth") && !parameters->value("snowSurfaceDampingDepth").toString().isEmpty())
{
snowModel.snowParameters.snowSurfaceDampingDepth = parameters->value("snowSurfaceDampingDepth").toDouble();
}
parameters->endGroup();
}
}
return true;
}


bool Crit3DProject::writeCriteria3DParameters()
{
QString fileName = getCompleteFileName(parametersFileName, PATH_SETTINGS);
Expand All @@ -538,6 +543,7 @@ bool Crit3DProject::writeCriteria3DParameters()
parameters->setValue("snow/skinThickness", snowModel.snowParameters.skinThickness);
parameters->setValue("snow/snowVegetationHeight", snowModel.snowParameters.snowVegetationHeight);
parameters->setValue("snow/soilAlbedo", snowModel.snowParameters.soilAlbedo);
parameters->setValue("snow/snowSurfaceDampingDepth", snowModel.snowParameters.snowSurfaceDampingDepth);

parameters->sync();
return true;
Expand Down
19 changes: 19 additions & 0 deletions bin/CRITERIA3D/dialogSnowSettings.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,13 @@ DialogSnowSettings::DialogSnowSettings(QWidget *parent) : QDialog(parent)
soilAlbedoValue->setValidator(doubleAlbedoVal);
soilAlbedoValue->setFixedWidth(70);

QLabel *dampingDepth = new QLabel(tr("Snow surface damping depth [m] "));
snowDampingDepthValue = new QLineEdit();
QDoubleValidator* doubleDampingDepthVal = new QDoubleValidator(0.0, 1.0, 2, snowDampingDepthValue);
doubleDampingDepthVal->setNotation(QDoubleValidator::StandardNotation);
snowDampingDepthValue->setValidator(doubleAlbedoVal);
snowDampingDepthValue->setFixedWidth(70);

layoutSettings->addWidget(rainfallThreshold, 0 , 0);
layoutSettings->addWidget(rainfallThresholdValue, 0 , 1);
layoutSettings->addWidget(snowThreshold, 1 , 0);
Expand All @@ -64,6 +71,8 @@ DialogSnowSettings::DialogSnowSettings(QWidget *parent) : QDialog(parent)
layoutSettings->addWidget(vegetationHeightValue, 4 , 1);
layoutSettings->addWidget(soilAlbedo, 5 , 0);
layoutSettings->addWidget(soilAlbedoValue, 5 , 1);
layoutSettings->addWidget(dampingDepth, 6 , 0);
layoutSettings->addWidget(snowDampingDepthValue, 6 , 1);

QDialogButtonBox *buttonBox = new QDialogButtonBox(QDialogButtonBox::Ok
| QDialogButtonBox::Cancel);
Expand Down Expand Up @@ -269,3 +278,13 @@ void DialogSnowSettings::setSoilAlbedoValue(double value)
{
soilAlbedoValue->setText(QLocale().toString(value));
}

double DialogSnowSettings::getSnowDampingDepthValue() const
{
return QLocale().toDouble(snowDampingDepthValue->text());
}

void DialogSnowSettings::setSnowDampingDepthValue(double value)
{
snowDampingDepthValue->setText(QLocale().toString(value));
}
4 changes: 4 additions & 0 deletions bin/CRITERIA3D/dialogSnowSettings.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,9 @@
double getSoilAlbedoValue() const;
void setSoilAlbedoValue(double value);

double getSnowDampingDepthValue() const;
void setSnowDampingDepthValue(double value);

bool checkEmptyValues();
bool checkWrongValues();
void accept();
Expand All @@ -39,6 +42,7 @@
QLineEdit *surfaceThickValue;
QLineEdit *vegetationHeightValue;
QLineEdit *soilAlbedoValue;
QLineEdit *snowDampingDepthValue;
};

#endif // DIALOGSNOWSETTINGS_H
24 changes: 11 additions & 13 deletions bin/CRITERIA3D/mainwindow.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2123,26 +2123,24 @@ void MainWindow::on_actionSnow_settings_triggered()
dialogSnowSetting.setSurfaceThickValue(myProject.snowModel.snowParameters.skinThickness);
dialogSnowSetting.setVegetationHeightValue(myProject.snowModel.snowParameters.snowVegetationHeight);
dialogSnowSetting.setSoilAlbedoValue(myProject.snowModel.snowParameters.soilAlbedo);
dialogSnowSetting.setSnowDampingDepthValue(myProject.snowModel.snowParameters.snowSurfaceDampingDepth);


dialogSnowSetting.exec();
if (dialogSnowSetting.result() != QDialog::Accepted)
{
return;
}
else
{
double tempMaxWithSnow = dialogSnowSetting.getRainfallThresholdValue();
double tempMinWithRain = dialogSnowSetting.getSnowThresholdValue();
double snowWaterHoldingCapacity = dialogSnowSetting.getWaterHoldingValue();
double skinThickness = dialogSnowSetting.getSurfaceThickValue();
double snowVegetationHeight = dialogSnowSetting.getVegetationHeightValue();
double soilAlbedo = dialogSnowSetting.getSoilAlbedoValue();
myProject.snowModel.snowParameters.tempMinWithRain = tempMinWithRain;
myProject.snowModel.snowParameters.tempMaxWithSnow = tempMaxWithSnow;
myProject.snowModel.snowParameters.snowWaterHoldingCapacity = snowWaterHoldingCapacity;
myProject.snowModel.snowParameters.skinThickness = skinThickness;
myProject.snowModel.snowParameters.snowVegetationHeight = snowVegetationHeight;
myProject.snowModel.snowParameters.soilAlbedo = soilAlbedo;
{
myProject.snowModel.snowParameters.tempMinWithRain = dialogSnowSetting.getSnowThresholdValue();
myProject.snowModel.snowParameters.tempMaxWithSnow = dialogSnowSetting.getRainfallThresholdValue();
myProject.snowModel.snowParameters.snowWaterHoldingCapacity = dialogSnowSetting.getWaterHoldingValue();
myProject.snowModel.snowParameters.skinThickness = dialogSnowSetting.getSurfaceThickValue();
myProject.snowModel.snowParameters.snowVegetationHeight = dialogSnowSetting.getVegetationHeightValue();
myProject.snowModel.snowParameters.soilAlbedo = dialogSnowSetting.getSoilAlbedoValue();
myProject.snowModel.snowParameters.snowSurfaceDampingDepth = dialogSnowSetting.getSnowDampingDepthValue();

if (!myProject.writeCriteria3DParameters())
{
myProject.logError("Error writing snow parameters");
Expand Down

0 comments on commit 32182c7

Please sign in to comment.