diff --git a/DATA/PROJECT/Troy/SETTINGS/parameters.ini b/DATA/PROJECT/Troy/SETTINGS/parameters.ini index 430bb9477..e6494ce56 100644 --- a/DATA/PROJECT/Troy/SETTINGS/parameters.ini +++ b/DATA/PROJECT/Troy/SETTINGS/parameters.ini @@ -73,3 +73,4 @@ snowSkinThickness=0.003 snowVegetationHeight=0.1 soilAlbedo=0.2 skinThickness=0.02 +snowSurfaceDampingDepth=0.05 diff --git a/agrolib/snow/snow.cpp b/agrolib/snow/snow.cpp index 5f6e5c497..171d1433b 100644 --- a/agrolib/snow/snow.cpp +++ b/agrolib/snow/snow.cpp @@ -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 } @@ -145,7 +145,6 @@ void Crit3DSnow::computeSnowBrooksModel() { double solarRadTot; double cloudCover; /*!< [-] */ - double prevIceContent, prevLWaterContent; /*!< [mm] */ double currentRatio; double airActualVapDensity; /*!< [kg m-3] */ @@ -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())) @@ -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; @@ -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) @@ -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); @@ -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; } diff --git a/agrolib/snow/snow.h b/agrolib/snow/snow.h index 8ec4183eb..0a8d0a6bc 100644 --- a/agrolib/snow/snow.h +++ b/agrolib/snow/snow.h @@ -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] */ @@ -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(); @@ -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 diff --git a/agrolib/snow/snowMaps.cpp b/agrolib/snow/snowMaps.cpp index cf9c17c31..4c79993d0 100644 --- a/agrolib/snow/snowMaps.cpp +++ b/agrolib/snow/snowMaps.cpp @@ -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.)); diff --git a/bin/CRITERIA3D/criteria3DProject.cpp b/bin/CRITERIA3D/criteria3DProject.cpp index b2e802e8f..f582086bd 100644 --- a/bin/CRITERIA3D/criteria3DProject.cpp +++ b/bin/CRITERIA3D/criteria3DProject.cpp @@ -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); @@ -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; diff --git a/bin/CRITERIA3D/dialogSnowSettings.cpp b/bin/CRITERIA3D/dialogSnowSettings.cpp index 5b69c0489..1800d6b61 100644 --- a/bin/CRITERIA3D/dialogSnowSettings.cpp +++ b/bin/CRITERIA3D/dialogSnowSettings.cpp @@ -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); @@ -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); @@ -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)); +} diff --git a/bin/CRITERIA3D/dialogSnowSettings.h b/bin/CRITERIA3D/dialogSnowSettings.h index 528b98c9d..3df5b6c5c 100644 --- a/bin/CRITERIA3D/dialogSnowSettings.h +++ b/bin/CRITERIA3D/dialogSnowSettings.h @@ -28,6 +28,9 @@ double getSoilAlbedoValue() const; void setSoilAlbedoValue(double value); + double getSnowDampingDepthValue() const; + void setSnowDampingDepthValue(double value); + bool checkEmptyValues(); bool checkWrongValues(); void accept(); @@ -39,6 +42,7 @@ QLineEdit *surfaceThickValue; QLineEdit *vegetationHeightValue; QLineEdit *soilAlbedoValue; + QLineEdit *snowDampingDepthValue; }; #endif // DIALOGSNOWSETTINGS_H diff --git a/bin/CRITERIA3D/mainwindow.cpp b/bin/CRITERIA3D/mainwindow.cpp index 7051131d1..88f2ebb21 100644 --- a/bin/CRITERIA3D/mainwindow.cpp +++ b/bin/CRITERIA3D/mainwindow.cpp @@ -2123,6 +2123,8 @@ 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) @@ -2130,19 +2132,15 @@ void MainWindow::on_actionSnow_settings_triggered() 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");