From a72e9fdfe266ca7fb57e87687eb6c34a48e3d886 Mon Sep 17 00:00:00 2001 From: Markus Schnorbus Date: Thu, 19 Sep 2024 17:04:57 -0700 Subject: [PATCH 1/9] Update .gitignore for local PCIC environment. --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) mode change 100644 => 100755 .gitignore diff --git a/.gitignore b/.gitignore old mode 100644 new mode 100755 index e0978408..39c479d1 --- a/.gitignore +++ b/.gitignore @@ -41,3 +41,6 @@ Makefile # Runtime artifacts Raven Raven_errors.txt + +# Directories +test/ From 0b3c2d2efb7a7ff0a406717128d719d91985b312 Mon Sep 17 00:00:00 2001 From: Markus Schnorbus Date: Thu, 19 Sep 2024 17:10:12 -0700 Subject: [PATCH 2/9] Initial changes to CEnthalpyModel::RouteMassinReservoir, Reservoir.cpp and Reservoir.h to introduce two lake layers to lake enthalpy model. --- src/EnergyTransport.cpp | 88 ++++++++++++++++++++++++++++++----------- src/Makefile | 8 ++-- src/Reservoir.cpp | 20 +++++++++- src/Reservoir.h | 5 +++ 4 files changed, 92 insertions(+), 29 deletions(-) mode change 100644 => 100755 src/EnergyTransport.cpp mode change 100644 => 100755 src/Makefile mode change 100644 => 100755 src/Reservoir.cpp mode change 100644 => 100755 src/Reservoir.h diff --git a/src/EnergyTransport.cpp b/src/EnergyTransport.cpp old mode 100644 new mode 100755 index 9c78c7d0..f34a0c04 --- a/src/EnergyTransport.cpp +++ b/src/EnergyTransport.cpp @@ -270,17 +270,21 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub double A_new = pRes->GetSurfaceArea(); double A_old = pRes->GetOldSurfaceArea(); double A_avg = 0.5 * (A_new + A_old); + double V_h_new = pRes->GetHypolimnionStorage (); + double V_h_old = pRes->GetOldHypolimnionStorage (); CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); double Acorr=1.0; double SW(0), LW(0), LW_in(0), temp_air(0), AET(0); + double hstar(0),ksed(0),Vsed=0.001; double T_new =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_new); double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p] / V_old); double Ts_new=ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / Vsed ); double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p] / Vsed ); + if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input @@ -300,7 +304,7 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub } Q_sens =hstar* A_avg * (temp_air -0.5*(T_new+T_old)); - Q_cond =ksed * A_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); + Q_cond =kdiff * A_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); Q_sw_in =(SW )*A_avg; Q_lw_in =(LW_in )*A_avg; Q_lw_out=(LW )*A_avg; @@ -340,15 +344,27 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // double tstep=Options.timestep; if((V_old<=0.0) || (V_new<=0.0)) { Res_mass=ResSedMass=0.0; return;} //handles dried out reservoir/lake - + + double V_h_new = pRes->GetHypolimnionStorage (); + double V_h_old = pRes->GetOldHypolimnionStorage (); + double V_e_new = V_new-V_h_new; + double V_e_old = V_old-V_h_old; + double Q_dn_new = 0; + double Q_dn_old = 0; + double Q_up_new = 0; + double Q_up_old = 0; + double A_h_new = pRes->GetMixingArea(); + double A_h_old = pRes->GetOldMixingArea(); + double A_h_avg = 0.5 * (A_h_new + A_h_old); + CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); double Acorr=1.0; - double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p]/V_old); + double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p]/V_e_old); double SW(0), LW(0), LW_in(0), temp_air(0), AET(0); - double hstar(0), ksed(0), Vsed=0.001; + double hstar(0); if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input Acorr =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area @@ -359,14 +375,38 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // LW =-STEFAN_BOLTZ*EMISS_WATER*pow(T_old+ZERO_CELSIUS,4); //[MJ/m2/d] //TMP DEBUG -time-lagged - should include in the N-R formulation AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] - - Vsed =pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; +// Vsed =pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; hstar =pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] - ksed =pRes->GetLakebedConductivity()/0.5/ pRes->GetLakebedThickness();//[MJ/m2/d/K] +// ksed =pRes->GetLakebedConductivity()/0.5/ pRes->GetLakebedThickness();//[MJ/m2/d/K] } - double T_sed_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p]/Vsed); - + double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p]/V_h_old); + + + double dens_e,dens_h,vdiff,kdiff; + + double a0= 999.842594; + double a1= 6.793952e-2; + double a2=-9.095290e-3; + double a3= 1.001685e-4; + double a4=-1.120083e-6; + double a5= 6.536332e-9; + + dens_e = a0+a1*T_old + a2*pow(T_old,2) + a3*pow(T_old,3) +a4*pow(T_old,4)+a5*pow(T_old,5); + dens_h = a0+a1*Ts_old + a2*pow(Ts_old,2) + a3*pow(Ts_old,3) +a4*pow(Ts_old,4)+a5*pow(Ts_old,5); + + double N2; + + if(dens_e tolerance) { - T_guess = ConvertVolumetricEnthalpyToTemperature(E_guess/V_new); - Ts_guess= ConvertVolumetricEnthalpyToTemperature(Es_guess/Vsed); + T_guess = ConvertVolumetricEnthalpyToTemperature(E_guess/V_e_new); + Ts_guess= ConvertVolumetricEnthalpyToTemperature(Es_guess/V_h_new); A[0][0] = 1.0 + 0.5 * tstep * Q_new/V_new; A[1][1] = 1.0; @@ -418,19 +458,19 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // A[0][1] = 0.0; if (E_guess!=0.0){ A[0][0]+= 0.5*tstep*hstar*A_new*T_guess/E_guess; - A[0][0]+= 0.5*tstep*ksed *A_avg*T_guess/E_guess; //[MJ/m2/d/K]*[m2]/[MJ/m3/K] = [m3/d] - A[1][0]-= 0.5*tstep*ksed *A_avg*T_guess/E_guess; + A[0][0]+= 0.5*tstep*kdiff *A_avg*T_guess/E_guess; //[MJ/m2/d/K]*[m2]/[MJ/m3/K] = [m3/d] + A[1][0]-= 0.5*tstep*kdiff *A_avg*T_guess/E_guess; } if (Es_guess!=0.0){ - A[0][1]-= 0.5*tstep*ksed *A_avg*Ts_guess/Es_guess; - A[1][1]+= 0.5*tstep*ksed *A_avg*Ts_guess/Es_guess; + A[0][1]-= 0.5*tstep*kdiff *A_avg*Ts_guess/Es_guess; + A[1][1]+= 0.5*tstep*kdiff *A_avg*Ts_guess/Es_guess; } //J_ij=A_ij+E*dA_ij/dE = Jacobian - J[0][0]=A[0][0]+0.5*tstep*(hstar*A_new+ksed*A_avg) * (TemperatureEnthalpyDerivative(E_guess /V_new)/ V_new); - J[0][1]=A[0][1]-0.5*tstep*( ksed*A_avg) * (TemperatureEnthalpyDerivative(Es_guess/Vsed )/ Vsed ); - J[1][0]=A[1][0]-0.5*tstep*( ksed*A_avg) * (TemperatureEnthalpyDerivative(E_guess /V_new)/ V_new); - J[1][1]=A[1][1]+0.5*tstep*( ksed*A_avg) * (TemperatureEnthalpyDerivative(Es_guess/Vsed )/ Vsed ); + J[0][0]=A[0][0]+0.5*tstep*(hstar*A_new+kdiff*A_avg) * (TemperatureEnthalpyDerivative(E_guess /V_e_new)/ V_e_new); + J[0][1]=A[0][1]-0.5*tstep*( kdiff*A_avg) * (TemperatureEnthalpyDerivative(Es_guess/V_h_new )/ V_h_new ); + J[1][0]=A[1][0]-0.5*tstep*( kdiff*A_avg) * (TemperatureEnthalpyDerivative(E_guess /V_e_new)/ V_e_new); + J[1][1]=A[1][1]+0.5*tstep*( kdiff*A_avg) * (TemperatureEnthalpyDerivative(Es_guess/V_h_new )/ V_h_new ); R[0] =-A[0][0]*E_guess-A[0][1]*Es_guess+B[0]; R[1] =-A[1][0]*E_guess-A[1][1]*Es_guess+B[1]; @@ -440,7 +480,7 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // dE = (J[0][1]*R[1]-J[1][1]*R[0])/den; dEs = (J[1][0]*R[0]-J[0][0]*R[1])/den; - change =sqrt(dE*dE+dEs*dEs)/ V_new / HCP_WATER; //convert to approx temp difference (for >0C water) + change =sqrt(dE*dE+dEs*dEs)/ V_e_new / HCP_WATER; //convert to approx temp difference (for >0C water) E_guess +=dE; Es_guess +=dEs; diff --git a/src/Makefile b/src/Makefile old mode 100644 new mode 100755 index 455c92a6..ee09ffe1 --- a/src/Makefile +++ b/src/Makefile @@ -25,12 +25,12 @@ LDFLAGS := CXXFLAGS += -std=c++11 -fPIC # OPTION 1) include netcdf - uncomment following two commands (assumes netCDF path = /usr/local): -#CXXFLAGS += -Dnetcdf -#LDLIBS += -L/usr/local -lnetcdf +CXXFLAGS += -Dnetcdf +LDLIBS += -L/usr/local -lnetcdf # OPTION 2) include lp_solve for water management optimization- uncomment following two commands (assumes liblpsolve55.a in ../lib/lp_solve_unix/ folder): -CXXFLAGS += -D_LPSOLVE_ -LDLIBS += -L../lib/lp_solve_unix -llpsolve55 +#CXXFLAGS += -D_LPSOLVE_ +#LDLIBS += -L../lib/lp_solve_unix -llpsolve55 # OPTION 3) if you use a OSX/BSD system, uncomment the LDFLAGS line below # this is to allow for use a 1Gb stack, see http://linuxtoosx.blogspot.ca/2010/10/stack-overflow-increasing-stack-limit.html diff --git a/src/Reservoir.cpp b/src/Reservoir.cpp old mode 100644 new mode 100755 index 3043dae3..32c5167f --- a/src/Reservoir.cpp +++ b/src/Reservoir.cpp @@ -18,7 +18,7 @@ void CReservoir::BaseConstructor(const string Name,const long SubID) _SBID=SubID; _pDownSB=NULL; - + // _mixing_depth=1; _lakebed_thick=1.0; _lakebed_cond =0.0; _lake_convcoeff=2.0; @@ -403,6 +403,14 @@ long CReservoir::GetSubbasinID () const { return _SBID; } // double CReservoir::GetStorage () const { return GetVolume(_stage); } +/// \returns reservoir hypoliminion storage [m3] +// +double CReservoir::GetHypolimnionStorage () const { return GetVolume(_stage-_mixing_depth); } + +/// \returns reservoir old hypoliminion storage [m3] +// +double CReservoir::GetOldHypolimnionStorage () const { return GetVolume(_stage_last-_mixing_depth); } + ////////////////////////////////////////////////////////////////// /// \returns current stage [m] // @@ -449,6 +457,16 @@ double CReservoir::GetSurfaceArea () const { return GetArea(_stage); } // double CReservoir::GetOldSurfaceArea () const { return GetArea(_stage_last); } +////////////////////////////////////////////////////////////////// +/// \returns current mixing area [m2] +// +double CReservoir::GetMixingArea () const { return GetArea(_stage-_mixing_depth); } + +////////////////////////////////////////////////////////////////// +/// \returns current old mixing area [m2] +// +double CReservoir::GetOldMixingArea () const { return GetArea(_stage_last-_mixing_depth); } + ////////////////////////////////////////////////////////////////// /// \returns lakebed thickness [m] // diff --git a/src/Reservoir.h b/src/Reservoir.h old mode 100644 new mode 100755 index 03894f57..36e8a7e8 --- a/src/Reservoir.h +++ b/src/Reservoir.h @@ -55,6 +55,7 @@ class CReservoir string _name; ///< reservoir name long _SBID; ///< subbasin ID double _max_capacity; ///< maximum reservoir capacity [m3] + double _mixing_depth; double _lakebed_thick; ///< lakebed thickness [m] double _lakebed_cond; ///< lakebed thermal conductivity [MJ/m/K/d] @@ -178,6 +179,8 @@ class CReservoir double GetStorage () const; //[m3] double GetOldStorage () const; //[m3] + double GetHypolimnionStorage () const; //[m3] + double GetOldHypolimnionStorage () const; //[m3] double GetOutflowRate () const; //[m3/s] double GetOldOutflowRate () const; //[m3/s] double GetIntegratedOutflow (const double &tstep) const; //[m3] @@ -191,6 +194,8 @@ class CReservoir double GetOldStage () const; //[m] double GetSurfaceArea () const; //[m2] double GetOldSurfaceArea () const; //[m2] + double GetMixingArea () const; //[m2] + double GetOldMixingArea () const; //[m2] double GetLakebedThickness () const; //[m] double GetLakebedConductivity () const; //[MJ/m/K/d] double GetLakeConvectionCoeff () const; //[MJ/m2/d/K] From 861a0702cd9d8b947fcc9fc15e890714784c8316 Mon Sep 17 00:00:00 2001 From: Markus Schnorbus Date: Fri, 20 Sep 2024 10:27:07 -0700 Subject: [PATCH 3/9] Additional modifications to CEnthalpyModel::RouteMassinReservoir to correct coefficient [A], constant [B], and Jacobian matrices. Also initialized value for _mixing_depth in CReservoir::BaseConstructor. --- src/EnergyTransport.cpp | 125 +++++++++++++++++++++------------------- src/Reservoir.cpp | 2 +- 2 files changed, 66 insertions(+), 61 deletions(-) diff --git a/src/EnergyTransport.cpp b/src/EnergyTransport.cpp index f34a0c04..8021d9f9 100755 --- a/src/EnergyTransport.cpp +++ b/src/EnergyTransport.cpp @@ -304,7 +304,7 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub } Q_sens =hstar* A_avg * (temp_air -0.5*(T_new+T_old)); - Q_cond =kdiff * A_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); + Q_cond =ksed * A_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); Q_sw_in =(SW )*A_avg; Q_lw_in =(LW_in )*A_avg; Q_lw_out=(LW )*A_avg; @@ -314,7 +314,7 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub return -(Q_sens + Q_cond + Q_sw_in + Q_lw_in + Q_lw_out + Q_lat + Q_rain) * tstep; //[MJ] } ////////////////////////////////////////////////////////////////// -/// \brief Calculates Res_mass (total reservoir energ) and ResSedMass (total sediment energy) in basin with reservoir at end of time step +/// \brief Calculates Res_mass (total reservoir energy) and ResSedMass (total sediment energy) in basin with reservoir at end of time step /// /// \param p [in] subbasin index /// \param aMout_new[] [in] Array of energy outflows at downstream end of each segment at end of current timestep [MJ/d] [size: nsegs] @@ -324,7 +324,7 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub /// \param tt [in] Time structure // void CEnthalpyModel::RouteMassInReservoir (const int p, // SB index - const double *aMout_new, // [MJ/d][size: nsegs ] + const double *aMout_new, // [MJ/d][size: nsegs] double &Res_mass, // [MJ] double &ResSedMass, // [MJ] const optStruct &Options, @@ -334,28 +334,30 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // CReservoir* pRes = _pModel->GetSubBasin(p)->GetReservoir(); int nSegments = _pModel->GetSubBasin(p)->GetNumSegments(); - double V_new=pRes->GetStorage (); - double V_old=pRes->GetOldStorage (); - double Q_new=pRes->GetOutflowRate ()*SEC_PER_DAY; - double Q_old=pRes->GetOldOutflowRate()*SEC_PER_DAY; - double A_new=pRes->GetSurfaceArea(); - double A_old=pRes->GetOldSurfaceArea(); - double A_avg=0.5*(A_new+A_old); - double tstep=Options.timestep; - + double V_new=pRes->GetStorage(); + double V_old=pRes->GetOldStorage(); + if((V_old<=0.0) || (V_new<=0.0)) { Res_mass=ResSedMass=0.0; return;} //handles dried out reservoir/lake - double V_h_new = pRes->GetHypolimnionStorage (); - double V_h_old = pRes->GetOldHypolimnionStorage (); + double V_h_new = pRes->GetHypolimnionStorage(); + double V_h_old = pRes->GetOldHypolimnionStorage(); double V_e_new = V_new-V_h_new; double V_e_old = V_old-V_h_old; - double Q_dn_new = 0; - double Q_dn_old = 0; - double Q_up_new = 0; - double Q_up_old = 0; + double Q_new=pRes->GetOutflowRate()*SEC_PER_DAY; + double Q_old=pRes->GetOldOutflowRate()*SEC_PER_DAY; + double A_new=pRes->GetSurfaceArea(); + double A_old=pRes->GetOldSurfaceArea(); + double A_avg=0.5*(A_new+A_old); double A_h_new = pRes->GetMixingArea(); double A_h_old = pRes->GetOldMixingArea(); double A_h_avg = 0.5 * (A_h_new + A_h_old); + + double Q_dn_new = 0.0; + double Q_dn_old = 0.0; + double Q_up_new = 0.0; + double Q_up_old = 0.0; + + double tstep=Options.timestep; CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); @@ -367,13 +369,10 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // double hstar(0); if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input Acorr =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area - temp_air =pHRU->GetForcingFunctions()->temp_ave; //[C] SW =pHRU->GetForcingFunctions()->SW_radia_net; //[MJ/m2/d] LW_in =pHRU->GetForcingFunctions()->LW_incoming; //[MJ/m2/d] - LW =-STEFAN_BOLTZ*EMISS_WATER*pow(T_old+ZERO_CELSIUS,4); //[MJ/m2/d] //TMP DEBUG -time-lagged - should include in the N-R formulation - AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] // Vsed =pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; hstar =pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] @@ -382,7 +381,7 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p]/V_h_old); - + // water density in each layer as a function of water temperature (Chapra, 2008) double dens_e,dens_h,vdiff,kdiff; double a0= 999.842594; @@ -395,23 +394,24 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // dens_e = a0+a1*T_old + a2*pow(T_old,2) + a3*pow(T_old,3) +a4*pow(T_old,4)+a5*pow(T_old,5); dens_h = a0+a1*Ts_old + a2*pow(Ts_old,2) + a3*pow(Ts_old,3) +a4*pow(Ts_old,4)+a5*pow(Ts_old,5); - double N2; - - if(dens_e= dens_h + if(dens_e0C water) + change =sqrt(dE*dE+dEs*dEs)/V_new/HCP_WATER; //convert to approx temp difference (for >0C water) E_guess +=dE; Es_guess +=dEs; @@ -494,6 +498,7 @@ double b = -3.1; delete [] R; for (int i=0;i<2;i++){delete [] A[i]; delete[] J[i]; } delete[] A; delete[] J; } + ////////////////////////////////////////////////////////////////// /// \brief returns total energy lost from subbasin reach over current time step [MJ] /// \param p [in] subbasin index diff --git a/src/Reservoir.cpp b/src/Reservoir.cpp index 32c5167f..abfb133e 100755 --- a/src/Reservoir.cpp +++ b/src/Reservoir.cpp @@ -18,7 +18,7 @@ void CReservoir::BaseConstructor(const string Name,const long SubID) _SBID=SubID; _pDownSB=NULL; - // _mixing_depth=1; + _mixing_depth=5.0; _lakebed_thick=1.0; _lakebed_cond =0.0; _lake_convcoeff=2.0; From 5fdcaa06aad5a30006ff8f193b8a4a50cee6596e Mon Sep 17 00:00:00 2001 From: Markus Schnorbus Date: Mon, 7 Oct 2024 10:33:46 -0700 Subject: [PATCH 4/9] Modified diffusive heat transfer between layers to estimate kdiff using procedure from Community Land Model v5. Added advective heat transfer (i.e. initializing Q_up_new and Q_dn_new). Some code clean-up. --- src/EnergyTransport.cpp | 280 +++++++++++++++++++++++++++------------- src/Reservoir.h | 2 +- src/Transport.h | 5 + 3 files changed, 193 insertions(+), 94 deletions(-) diff --git a/src/EnergyTransport.cpp b/src/EnergyTransport.cpp index 8021d9f9..0dc25513 100755 --- a/src/EnergyTransport.cpp +++ b/src/EnergyTransport.cpp @@ -270,8 +270,13 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub double A_new = pRes->GetSurfaceArea(); double A_old = pRes->GetOldSurfaceArea(); double A_avg = 0.5 * (A_new + A_old); - double V_h_new = pRes->GetHypolimnionStorage (); - double V_h_old = pRes->GetOldHypolimnionStorage (); + double V_h_new = pRes->GetHypolimnionStorage(); + double V_h_old = pRes->GetOldHypolimnionStorage(); + double V_e_new = V_new-V_h_new; + double V_e_old = V_old-V_h_old; + double A_h_new = pRes->GetMixingArea(); + double A_h_old = pRes->GetOldMixingArea(); + double A_h_avg = 0.5 * (A_h_new + A_h_old); CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); @@ -279,37 +284,38 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub double SW(0), LW(0), LW_in(0), temp_air(0), AET(0); - double hstar(0),ksed(0),Vsed=0.001; - double T_new =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_new); - double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p] / V_old); - double Ts_new=ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / Vsed ); - double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p] / Vsed ); + double hstar(0),ksed(0); //Vsed=0.001; + double T_new =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_e_new); + double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p] / V_e_old); + double Ts_new=ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_h_new ); + double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p] / V_h_old ); + if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input Acorr =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area - temp_air =pHRU->GetForcingFunctions()->temp_ave; //[C] - SW =pHRU->GetForcingFunctions()->SW_radia_net; //[MJ/m2/d] - not using canopy correction! - LW_in =pHRU->GetForcingFunctions()->LW_incoming; //[MJ/m2/d] + temp_air =pHRU->GetForcingFunctions()->temp_ave; //[C] + SW =pHRU->GetForcingFunctions()->SW_radia_net; //[MJ/m2/d] - not using canopy correction! + LW_in =pHRU->GetForcingFunctions()->LW_incoming; //[MJ/m2/d] LW =-STEFAN_BOLTZ*EMISS_WATER*0.5*(pow(T_new+ZERO_CELSIUS,4)+pow(T_old+ZERO_CELSIUS,4)); - AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] //*pHRU->GetArea()/A_avg + AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] //*pHRU->GetArea()/A_avg - hstar = pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] - Vsed = pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; - ksed = pRes->GetLakebedConductivity() / 0.5 / pRes->GetLakebedThickness();//[MJ/m2/d/K] + hstar = pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] + // Vsed = pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; + ksed = pRes->GetLakebedConductivity() / 0.5 / pRes->GetLakebedThickness(); //[MJ/m2/d/K] } - Q_sens =hstar* A_avg * (temp_air -0.5*(T_new+T_old)); - Q_cond =ksed * A_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); - Q_sw_in =(SW )*A_avg; - Q_lw_in =(LW_in )*A_avg; - Q_lw_out=(LW )*A_avg; - Q_lat =-(AET * DENSITY_WATER * LH_VAPOR)*A_avg; - Q_rain =_aMresRain[p]; + Q_sens =hstar* A_avg * (temp_air -0.5*(T_new+T_old)); + Q_cond =ksed * A_h_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); + Q_sw_in =(SW )*A_avg; + Q_lw_in =(LW_in )*A_avg; + Q_lw_out=(LW )*A_avg; + Q_lat =-(AET * DENSITY_WATER * LH_VAPOR)*A_avg; + Q_rain =_aMresRain[p]; return -(Q_sens + Q_cond + Q_sw_in + Q_lw_in + Q_lw_out + Q_lat + Q_rain) * tstep; //[MJ] } @@ -334,8 +340,8 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // CReservoir* pRes = _pModel->GetSubBasin(p)->GetReservoir(); int nSegments = _pModel->GetSubBasin(p)->GetNumSegments(); - double V_new=pRes->GetStorage(); - double V_old=pRes->GetOldStorage(); + double V_new = pRes->GetStorage(); + double V_old = pRes->GetOldStorage(); if((V_old<=0.0) || (V_new<=0.0)) { Res_mass=ResSedMass=0.0; return;} //handles dried out reservoir/lake @@ -343,69 +349,157 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // double V_h_old = pRes->GetOldHypolimnionStorage(); double V_e_new = V_new-V_h_new; double V_e_old = V_old-V_h_old; - double Q_new=pRes->GetOutflowRate()*SEC_PER_DAY; - double Q_old=pRes->GetOldOutflowRate()*SEC_PER_DAY; - double A_new=pRes->GetSurfaceArea(); - double A_old=pRes->GetOldSurfaceArea(); - double A_avg=0.5*(A_new+A_old); + double Q_new = pRes->GetOutflowRate() * SEC_PER_DAY; + double Q_old = pRes->GetOldOutflowRate() * SEC_PER_DAY; + double A_new = pRes->GetSurfaceArea(); + double A_old = pRes->GetOldSurfaceArea(); + double A_avg = 0.5 * (A_new + A_old); double A_h_new = pRes->GetMixingArea(); double A_h_old = pRes->GetOldMixingArea(); double A_h_avg = 0.5 * (A_h_new + A_h_old); - double Q_dn_new = 0.0; + double tstep=Options.timestep; + + + double Q_dn_new(0),Q_up_new(0); double Q_dn_old = 0.0; - double Q_up_new = 0.0; double Q_up_old = 0.0; - double tstep=Options.timestep; - - CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); + double Q_vert = (V_h_new - V_h_old)/tstep; + if (Q_vert > 0){ + Q_dn_new = Q_vert; Q_up_new = 0.0; + } + else + { + Q_dn_new = 0.0; Q_up_new = -Q_vert; + } + + CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); double Acorr=1.0; - double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p]/V_e_old); + double Ts_old =ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p]/V_h_old); double SW(0), LW(0), LW_in(0), temp_air(0), AET(0); double hstar(0); +//Vsed=0.001; + double u2(0); + if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input - Acorr =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area - temp_air =pHRU->GetForcingFunctions()->temp_ave; //[C] - SW =pHRU->GetForcingFunctions()->SW_radia_net; //[MJ/m2/d] - LW_in =pHRU->GetForcingFunctions()->LW_incoming; //[MJ/m2/d] - LW =-STEFAN_BOLTZ*EMISS_WATER*pow(T_old+ZERO_CELSIUS,4); //[MJ/m2/d] //TMP DEBUG -time-lagged - should include in the N-R formulation - AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] + Acorr =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area + temp_air =pHRU->GetForcingFunctions()->temp_ave; //[C] + SW =pHRU->GetForcingFunctions()->SW_radia_net; //[MJ/m2/d] + LW_in =pHRU->GetForcingFunctions()->LW_incoming; //[MJ/m2/d] + LW =-STEFAN_BOLTZ*EMISS_WATER*pow(T_old+ZERO_CELSIUS,4); //[MJ/m2/d] //TMP DEBUG -time-lagged - should include in the N-R formulation + AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] // Vsed =pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; - hstar =pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] + hstar =pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] // ksed =pRes->GetLakebedConductivity()/0.5/ pRes->GetLakebedThickness();//[MJ/m2/d/K] } - double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p]/V_h_old); + + +///////////// Option 1 to calculate kdiff +//// water density in each layer as a function of water temperature (Chapra, 2008) +//double dens_e,dens_h,vdiff,kdiff; +//double a0= 999.842594; +//double a1= 6.793952e-2; +//double a2=-9.095290e-3; +//double a3= 1.001685e-4; +//double a4=-1.120083e-6; +//double a5= 6.536332e-9; +//dens_e = a0+a1*T_old + a2*pow(T_old,2) + a3*pow(T_old,3) +a4*pow(T_old,4)+a5*pow(T_old,5); +//dens_h = a0+a1*Ts_old + a2*pow(Ts_old,2) + a3*pow(Ts_old,3) +a4*pow(Ts_old,4)+a5*pow(Ts_old,5); +//// Brunt-Väisälä frequency at thermocline depth (assume dz = 1.0m) (Jennings et al. 2012) +//double N2=1.0e-12; // value when dens_e >= dens_h +//if(dens_e kdiff = md*(Ke+Ked+Km); + u2 = pHRU->GetForcingFunctions()->wind_vel; - // water density in each layer as a function of water temperature (Chapra, 2008) - double dens_e,dens_h,vdiff,kdiff; + double kstar = 0.0; + double wstar = 0.0; //Surface friction velocity [m/s] + double Ri = 0.0; //Richardson number + double N2 = 5.0E-12; //Brunt-Väisälä frequency when dens_e >= dens_h [1/s2] + double dens_e, dens_h; //Water density in each layer [kg/m3] + double Ke(0.0),Ked(0.0),Km(0.0),kdiff(0.0); //Diffusion coefficient [MJ/m2/K/d] + double md = 1.0; //Diffusion scaling factor in each layer + + //Water density in each layer as a function of water temperature. Reference: (Hostetler and Bartlein, 1990) + double a0 = -1.954e-5; + dens_e = DENSITY_WATER*(1.0+a0*pow(abs((T_old + ZERO_CELSIUS)-277.0),1.68)); + dens_h = DENSITY_WATER*(1.0+a0*pow(abs((Ts_old + ZERO_CELSIUS)-277.0),1.68)); - double a0= 999.842594; - double a1= 6.793952e-2; - double a2=-9.095290e-3; - double a3= 1.001685e-4; - double a4=-1.120083e-6; - double a5= 6.536332e-9; + //Estimate layer thicknesses + double Ze = V_e_old/A_old; //Average epilimnion thickness [m] + double Zh = V_h_old/A_h_old; //Average hppolimnion thickness [m] + double Zc = (Ze+Zh)/2.0; //Characteristic length [m] - dens_e = a0+a1*T_old + a2*pow(T_old,2) + a3*pow(T_old,3) +a4*pow(T_old,4)+a5*pow(T_old,5); - dens_h = a0+a1*Ts_old + a2*pow(Ts_old,2) + a3*pow(Ts_old,3) +a4*pow(Ts_old,4)+a5*pow(Ts_old,5); + //Molecular diffusion coefficient, Km can be calculated as per Community Land Model v5 (Lawrence et al, 2020) + Km = TC_WATER/HCP_WATER; //[MJ/m2/K/d] - // Brunt-Väisälä frequency at thermocline depth (assume dz = 1.0m) (Jennings et al. 2012) - double N2=1.0e-12; // value when dens_e >= dens_h + //Estimate density gradient dp/dz over length Zc; ===> Brunt-Väisälä frequency across thermocline depth (Jennings et al. 2012) if(dens_e FREEZING_TEMP){ + wstar = 0.0012*u2; //Surfaec friction velocity [m/s] + kstar = 6.6*pow(u2, -1.84) * sqrt(abs(sin(pHRU->GetLatRad()))); + Ri = 40.0*N2*pow(VON_KARMAN,2.0)*pow(Zc,2.0); + Ri /= pow(wstar,2.0); + Ri /= exp(-2*VON_KARMAN*Zc); + Ri += 1.0; + Ri = (-1 + sqrt(Ri))/20.0; + Ke = VON_KARMAN*wstar*Zc/(1.0 + 37.0*pow(Ri,2.0)) * exp(-VON_KARMAN*Zc); //[MJ/m2/K/s] + Ke *= SEC_PER_DAY; //[MJ/m2/K/d] + } + + //Enhanced diffusivity intended to represent unresolved mixing processes; ====> Community Land Model v5 (Lawrence et al, 2020) + if (N2 >= 7.5e-5){ + Ked = 1.04e-8 * pow(N2,-0.43)*SEC_PER_DAY; //[MJ/m2/K/d] + } + //Increase the overall diffusivity for large lakes, intended to represent + //3-dimensional mixing processes such as cased by horizontal temperature; Gradient can be calculated as per Community Land Model v5 (Lawrence et al, 2020) + if ((Ze+Zh) >= 25.0){ + md = 10.0; //Diffusion scaling factor in each layer + } + + //Final diffusion coefficient + kdiff = md*(Ke+Ked+Km)*HCP_WATER; //[MJ/m2/K/d] + + + //TEMPORARY OUPUT ====> To test which parameter is abnormal + //cout<<"Ze = "<0C water) + + change =sqrt(dE*dE+dEs*dEs)/V_new/HCP_WATER; //convert to approx temp difference (for >0C water) E_guess +=dE; Es_guess +=dEs; @@ -1218,8 +1310,10 @@ void CEnthalpyModel::WriteMinorOutput(const optStruct& Options,const time_struct if ((pRes != NULL) && (pRes->GetHRUIndex() != DOESNT_EXIST)) { double HRUarea= _pModel->GetHydroUnit(pRes->GetHRUIndex())->GetArea()*M2_PER_KM2; - double V_sed = pRes->GetLakebedThickness() * HRUarea; + //double V_sed = pRes->GetLakebedThickness() * HRUarea; double V_new = pRes->GetStorage(); + double V_h_new = pRes->GetHypolimnionStorage(); + double V_e_new = V_new-V_h_new; mult = 1.0 / HRUarea; @@ -1230,20 +1324,20 @@ void CEnthalpyModel::WriteMinorOutput(const optStruct& Options,const time_struct GetEnergyLossesFromLake(p,Q_sens,Q_cond,Q_lat,Q_sw_in,Q_lw_in,Q_lw_out,Q_rain); - double lakeTemp =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_new); - double sedTemp =ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_sed); + double lakeTemp =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_e_new); + double sedTemp =ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_h_new); double inflowTemp=ConvertVolumetricEnthalpyToTemperature( Ein / Qin ); - double pctFroz =ConvertVolumetricEnthalpyToIceContent (_aMres[p] / V_new); + double pctFroz =ConvertVolumetricEnthalpyToIceContent (_aMres[p] / V_e_new); _LAKEOUT << mult * Ein << "," << mult * Eout << ","; _LAKEOUT << mult * Q_rain << "," << mult * Q_sens << ","; _LAKEOUT << mult * Q_cond << "," << mult * Q_lat << ","; _LAKEOUT << mult * Q_sw_in << "," << mult * Q_lw_in << "," << mult * Q_lw_out << ","; _LAKEOUT << mult * _aMres[p] << ","; - if (V_new!=0){_LAKEOUT << lakeTemp << ",";}else{_LAKEOUT<<",";} - if (V_sed!=0){_LAKEOUT << sedTemp << ",";}else{_LAKEOUT<<",";} - if (Qin !=0){_LAKEOUT << inflowTemp << ",";}else{_LAKEOUT<<",";} - if (V_new!=0){_LAKEOUT << pctFroz << ",";}else{_LAKEOUT<<",";} + if (V_e_new!=0){_LAKEOUT << lakeTemp << ",";}else{_LAKEOUT<<",";} + if (V_h_new!=0){_LAKEOUT << sedTemp << ",";}else{_LAKEOUT<<",";} + if (Qin !=0){_LAKEOUT << inflowTemp << ",";}else{_LAKEOUT<<",";} + if (V_e_new!=0){_LAKEOUT << pctFroz << ",";}else{_LAKEOUT<<",";} } } } diff --git a/src/Reservoir.h b/src/Reservoir.h index 36e8a7e8..45fb8528 100755 --- a/src/Reservoir.h +++ b/src/Reservoir.h @@ -195,7 +195,7 @@ class CReservoir double GetSurfaceArea () const; //[m2] double GetOldSurfaceArea () const; //[m2] double GetMixingArea () const; //[m2] - double GetOldMixingArea () const; //[m2] + double GetOldMixingArea () const; //[m2] double GetLakebedThickness () const; //[m] double GetLakebedConductivity () const; //[MJ/m/K/d] double GetLakeConvectionCoeff () const; //[MJ/m2/d/K] diff --git a/src/Transport.h b/src/Transport.h index 4d0514a0..4b7c0178 100644 --- a/src/Transport.h +++ b/src/Transport.h @@ -217,6 +217,11 @@ class CConstituentModel double *_aMout_res; ///< array storing reservoir mass outflow [mg/d] or heat loss [MJ/d] [size: nSubBasins] double *_aMout_res_last; ///< array storing reservoir mass outflow [mg/d] or enthalpy outflow [MJ/d] at start of timestep [size: nSubBasins] double *_aMresRain; ///< array storing reservoir rain inputs [mg/d] or enthalpy input [MJ/d] [size: nSubBasins] + + double *_aMres_epo; ///< array used for storing epolimnion layer masses [mg] or enthalpy [MJ] [size: nSubBasins] + double *_aMres_epo_last; ///< array used for storing epolimnion layer masses [mg] at start of timestep [size: nSubBasins] + double *_aMres_hyp; ///< array used for storing hypolimnion layer masses [mg] or enthalpy [MJ] [size: nSubBasins] + double *_aMres_hyp_last; ///< array used for storing hypolimnion layer masses [mg] at start of timestep [size: nSubBasins] // Mass balance tracking variables double *_channel_storage; ///< array storing channel storage [mg] or [MJ] [size: nSubBasins] From edce19bc6b75e8795153a5300026dcbbb9e09e7e Mon Sep 17 00:00:00 2001 From: Markus Schnorbus Date: Thu, 10 Oct 2024 10:31:47 -0700 Subject: [PATCH 5/9] Added vertical heat advection between lake layers in CENthalpyModel::RouteMassinReservoir --- src/EnergyTransport.cpp | 28 +++++++++++++++------------- src/Reservoir.cpp | 2 ++ src/Reservoir.h | 2 ++ 3 files changed, 19 insertions(+), 13 deletions(-) diff --git a/src/EnergyTransport.cpp b/src/EnergyTransport.cpp index 0dc25513..eb7a783c 100755 --- a/src/EnergyTransport.cpp +++ b/src/EnergyTransport.cpp @@ -361,19 +361,15 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // double tstep=Options.timestep; - double Q_dn_new(0),Q_up_new(0); - double Q_dn_old = 0.0; - double Q_up_old = 0.0; - - double Q_vert = (V_h_new - V_h_old)/tstep; - if (Q_vert > 0){ - Q_dn_new = Q_vert; Q_up_new = 0.0; - } - else - { - Q_dn_new = 0.0; Q_up_new = -Q_vert; - } - +// double Q_dn_new(0),Q_up_new(0),Q_dn_old(0),Q_up_old(0),Q_vert; + double Q_dn_new(0),Q_dn_old(0),Q_up_new(0),Q_up_old(0); + + double Q_vert = (V_h_new - V_h_old)/tstep; + if (Q_vert> 0){ + Q_dn_new = Q_vert; Q_up_new = 0.0; + }else{ + Q_dn_new = 0; Q_up_new = -Q_vert; + } CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); double Acorr=1.0; @@ -585,6 +581,12 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // } Res_mass =E_guess; ResSedMass=Es_guess; + + double _Q_dn_old = Q_dn_new; + double _Q_up_old = Q_up_new; + + cout<<"_Q_dn_old = "<<_Q_dn_old<<"\n"; + cout<<"_Q_up_old = "<<_Q_up_old<<"\n"; delete [] B; delete [] R; diff --git a/src/Reservoir.cpp b/src/Reservoir.cpp index abfb133e..08bdf0eb 100755 --- a/src/Reservoir.cpp +++ b/src/Reservoir.cpp @@ -18,6 +18,8 @@ void CReservoir::BaseConstructor(const string Name,const long SubID) _SBID=SubID; _pDownSB=NULL; + _Q_dn_old = 0.0; + _Q_up_old = 0.0; _mixing_depth=5.0; _lakebed_thick=1.0; _lakebed_cond =0.0; diff --git a/src/Reservoir.h b/src/Reservoir.h index 45fb8528..8ae38567 100755 --- a/src/Reservoir.h +++ b/src/Reservoir.h @@ -108,6 +108,8 @@ class CReservoir double _stage_last; ///< stage at beginning of current time step [m] double _Qout; ///< outflow corresponding to current stage [m3/s] double _Qout_last; ///< outflow at beginning of current time step [m3/s] + double _Q_dn_old; + double _Q_up_old; double _MB_losses; ///< losses over current time step [m3] double _AET; ///< losses through AET only [m3] double _Precip; ///< gains through precipitation [m3] From e1e1024575629780c9884cd476e48acc8814e545 Mon Sep 17 00:00:00 2001 From: Markus Schnorbus Date: Sun, 20 Oct 2024 12:27:58 -0700 Subject: [PATCH 6/9] Modifications and upgrades to the calculation of kdiff, including correcting error in calculation of Richardson number (replaced von Karman with k*) and adding upper limit on value of kdiff. Also added a function to get the mixing depth of the lake/reservoir. Some formatting changes to the code. --- src/EnergyTransport.cpp | 290 ++++++++++++++++++++++++---------------- src/EnergyTransport.h | 2 +- src/RavenInclude.h | 2 +- src/Reservoir.cpp | 7 +- src/Reservoir.h | 1 + 5 files changed, 183 insertions(+), 119 deletions(-) diff --git a/src/EnergyTransport.cpp b/src/EnergyTransport.cpp index eb7a783c..35a750bc 100755 --- a/src/EnergyTransport.cpp +++ b/src/EnergyTransport.cpp @@ -255,13 +255,14 @@ void CEnthalpyModel::SetHyporheicLayer(const int m) /// \param Q_rain [out] energy gain from precip inputs [MJ/d] /// \returns total energy lost from reach over current time step [MJ] // -double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, double &Q_cond, double &Q_lat, double &Q_sw_in, double &Q_lw_in, double &Q_lw_out, double &Q_rain) const +double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, double &Q_conv, double &Q_lat, double &Q_sw_in, double &Q_lw_in, double &Q_lw_out, double &Q_rain, double &Q_adv, double &Ri, double &N2, double &kdiff) const { double tstep = _pModel->GetOptStruct()->timestep; CSubBasin *pBasin = _pModel->GetSubBasin(p); CReservoir *pRes = _pModel->GetSubBasin(p)->GetReservoir(); if (pRes==NULL){return 0.0;} + CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); double V_new = pRes->GetStorage(); double V_old = pRes->GetOldStorage(); @@ -277,48 +278,117 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub double A_h_new = pRes->GetMixingArea(); double A_h_old = pRes->GetOldMixingArea(); double A_h_avg = 0.5 * (A_h_new + A_h_old); - - CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); + double zm = pRes->GetMixingDepth(); + + double Q_dn_new(0),Q_dn_old(0),Q_up_new(0),Q_up_old(0); + double Q_vert = (V_h_new - V_h_old)/tstep; + if (Q_vert> 0){ + Q_dn_new = Q_vert; Q_up_new = 0.0; + } else { + Q_dn_new = 0.0; Q_up_new = -Q_vert; + } double Acorr=1.0; - double SW(0), LW(0), LW_in(0), temp_air(0), AET(0); - - double hstar(0),ksed(0); //Vsed=0.001; + double hstar(0), ksed(0), Vsed=0.001; + double u2(0); + double T_new =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_e_new); double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p] / V_e_old); double Ts_new=ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_h_new ); double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p] / V_h_old ); - - if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input - Acorr =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area - temp_air =pHRU->GetForcingFunctions()->temp_ave; //[C] SW =pHRU->GetForcingFunctions()->SW_radia_net; //[MJ/m2/d] - not using canopy correction! LW_in =pHRU->GetForcingFunctions()->LW_incoming; //[MJ/m2/d] - LW =-STEFAN_BOLTZ*EMISS_WATER*0.5*(pow(T_new+ZERO_CELSIUS,4)+pow(T_old+ZERO_CELSIUS,4)); - AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] //*pHRU->GetArea()/A_avg - hstar = pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] - // Vsed = pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; + Vsed = pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; ksed = pRes->GetLakebedConductivity() / 0.5 / pRes->GetLakebedThickness(); //[MJ/m2/d/K] } - Q_sens =hstar* A_avg * (temp_air -0.5*(T_new+T_old)); - Q_cond =ksed * A_h_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); - Q_sw_in =(SW )*A_avg; - Q_lw_in =(LW_in )*A_avg; - Q_lw_out=(LW )*A_avg; - Q_lat =-(AET * DENSITY_WATER * LH_VAPOR)*A_avg; - Q_rain =_aMresRain[p]; + ///////////// Option 2 to calculate kdiff===> kdiff = md*(Ke+Ked+Km); + u2 = pHRU->GetForcingFunctions()->wind_vel; + + double kstar = 0.0; + double wstar = 0.0; //Surface friction velocity [m/s] + Ri = 0.0; //Richardson number + N2 = 5.0E-12; //Brunt-Väisälä frequency when dens_e >= dens_h [1/s2] + double dens_e, dens_h; //Water density in each layer [kg/m3] + double Ke(0.0),Ked(0.0),Km(0.0); //kdiff(0.0); //Diffusion coefficient [MJ/m2/K/d] + kdiff = 0.0; + double md = 1.0; //Diffusion scaling factor in each layer + + //Water density in each layer as a function of water temperature. Reference: + // (Hostetler and Bartlein, 1990) + double a0 = -1.954e-5; + dens_e = DENSITY_WATER*(1.0+a0*pow(abs((T_old + ZERO_CELSIUS)-277.0),1.68)); + dens_h = DENSITY_WATER*(1.0+a0*pow(abs((Ts_old + ZERO_CELSIUS)-277.0),1.68)); + + //Estimate layer thicknesses + double Ze = V_e_old/A_old; //Average epilimnion thickness [m] + double Zh = V_h_old/A_h_old; //Average hppolimnion thickness [m] + double Zc = (Ze+Zh)/2.0; //Characteristic length [m] + + //Molecular diffusion coefficient, Km can be calculated as per Community Land Model v5 + //(Lawrence et al, 2020) + Km = TC_WATER/HCP_WATER; //[MJ/m2/K/d] + + //Estimate density gradient dp/dz over length Zc; ===> Brunt-Väisälä frequency across + //the thermocline depth (Jennings et al. 2012) + if(dens_e FREEZING_TEMP){ + wstar = 0.0012*u2; //Surface friction velocity [m/s] + kstar = 6.6*pow(u2, -1.84) * sqrt(abs(sin(pHRU->GetLatRad()))); + Ri = 40.0*N2*pow(VON_KARMAN,2.0)*pow(zm,2.0); + Ri /= pow(wstar,2.0); + Ri /= exp(-2*kstar*zm); + Ri += 1.0; + Ri = (-1 + sqrt(Ri))/20.0; // Richardson number + if(Ri < 1e+20){ + Ke = VON_KARMAN*wstar*zm/(1.0 + 37.0*pow(Ri,2.0)) * exp(-kstar*zm); //[MJ/m2/K/s] + Ke *= SEC_PER_DAY; //[MJ/m2/K/d] + } + } + + //Enhanced diffusivity intended to represent unresolved mixing processes + if (N2 >= 7.5e-5){ + Ked = 1.04e-8 * pow(N2,-0.43)*SEC_PER_DAY; //[MJ/m2/K/d] + } + + //Increase the overall diffusivity for large lakes, intended to represent + //3-dimensional mixing processes such as cased by horizontal temperature; + // Gradient can be calculated as per Community Land Model v5 (Lawrence et al, 2020) + if ((Ze+Zh) >= 25.0){ + md = 10.0; //Diffusion scaling factor in each layer + } + + //Final diffusion coefficient + //kdiff = md*(Ke+Ked+Km)*HCP_WATER; //[MJ/m2/K/d] + kdiff = threshMin(md*(Ke+Ked+Km)*HCP_WATER, 5.0, 0.0); //[MJ/m2/K/d] + + + Q_sens =hstar * A_avg * (temp_air -0.5*(T_new+T_old)); + Q_conv =kdiff * A_h_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); + Q_sw_in =(SW )*A_avg; + Q_lw_in =(LW_in )*A_avg; + Q_lw_out=(LW )*A_avg; + Q_lat =-(AET * DENSITY_WATER * LH_VAPOR)*A_avg; + Q_rain =_aMresRain[p]; + Q_adv = (Q_up_new*_aMsed[p]/V_e_new - Q_dn_new*_aMres[p]/V_h_new) * tstep; - return -(Q_sens + Q_cond + Q_sw_in + Q_lw_in + Q_lw_out + Q_lat + Q_rain) * tstep; //[MJ] + //return -(Q_sens + Q_cond + Q_sw_in + Q_lw_in + Q_lw_out + Q_lat + Q_rain) * tstep; //[MJ] + return -(Q_sens + Q_sw_in + Q_lw_in + Q_lw_out + Q_lat + Q_rain) * tstep; //[MJ] } + ////////////////////////////////////////////////////////////////// /// \brief Calculates Res_mass (total reservoir energy) and ResSedMass (total sediment energy) in basin with reservoir at end of time step /// @@ -357,19 +427,17 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // double A_h_new = pRes->GetMixingArea(); double A_h_old = pRes->GetOldMixingArea(); double A_h_avg = 0.5 * (A_h_new + A_h_old); + double zm = pRes->GetMixingDepth(); double tstep=Options.timestep; - -// double Q_dn_new(0),Q_up_new(0),Q_dn_old(0),Q_up_old(0),Q_vert; double Q_dn_new(0),Q_dn_old(0),Q_up_new(0),Q_up_old(0); - - double Q_vert = (V_h_new - V_h_old)/tstep; + double Q_vert = (V_h_new - V_h_old)/tstep; if (Q_vert> 0){ - Q_dn_new = Q_vert; Q_up_new = 0.0; - }else{ - Q_dn_new = 0; Q_up_new = -Q_vert; - } + Q_dn_new = Q_vert; Q_up_new = 0.0; + } else { + Q_dn_new = 0; Q_up_new = -Q_vert; + } CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); double Acorr=1.0; @@ -377,8 +445,7 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // double Ts_old =ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p]/V_h_old); double SW(0), LW(0), LW_in(0), temp_air(0), AET(0); - double hstar(0); -//Vsed=0.001; + double hstar(0); // ksed=0.0, Vsed=0.001; double u2(0); if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input @@ -393,8 +460,6 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // // ksed =pRes->GetLakebedConductivity()/0.5/ pRes->GetLakebedThickness();//[MJ/m2/d/K] } - - ///////////// Option 1 to calculate kdiff //// water density in each layer as a function of water temperature (Chapra, 2008) //double dens_e,dens_h,vdiff,kdiff; @@ -404,8 +469,8 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // //double a3= 1.001685e-4; //double a4=-1.120083e-6; //double a5= 6.536332e-9; -//dens_e = a0+a1*T_old + a2*pow(T_old,2) + a3*pow(T_old,3) +a4*pow(T_old,4)+a5*pow(T_old,5); -//dens_h = a0+a1*Ts_old + a2*pow(Ts_old,2) + a3*pow(Ts_old,3) +a4*pow(Ts_old,4)+a5*pow(Ts_old,5); +//dens_e = a0+a1*T_old + a2*pow(T_old,2) + a3*pow(T_old,3) + a4*pow(T_old,4) + a5*pow(T_old,5); +//dens_h = a0+a1*Ts_old + a2*pow(Ts_old,2) + a3*pow(Ts_old,3) + a4*pow(Ts_old,4) + a5*pow(Ts_old,5); //// Brunt-Väisälä frequency at thermocline depth (assume dz = 1.0m) (Jennings et al. 2012) //double N2=1.0e-12; // value when dens_e >= dens_h //if(dens_e kdiff = md*(Ke+Ked+Km); + ///////////// Option 2 to calculate kdiff===> kdiff = md*(Ke+Ked+Km); u2 = pHRU->GetForcingFunctions()->wind_vel; double kstar = 0.0; - double wstar = 0.0; //Surface friction velocity [m/s] - double Ri = 0.0; //Richardson number - double N2 = 5.0E-12; //Brunt-Väisälä frequency when dens_e >= dens_h [1/s2] - double dens_e, dens_h; //Water density in each layer [kg/m3] - double Ke(0.0),Ked(0.0),Km(0.0),kdiff(0.0); //Diffusion coefficient [MJ/m2/K/d] - double md = 1.0; //Diffusion scaling factor in each layer - - //Water density in each layer as a function of water temperature. Reference: (Hostetler and Bartlein, 1990) + double wstar = 0.0; //Surface friction velocity [m/s] + double Ri = 0.0; //Richardson number + double N2 = 5.0E-12; //Brunt-Väisälä frequency when dens_e >= dens_h [1/s2] + double dens_e, dens_h; //Water density in each layer [kg/m3] + double Ke(0.0),Ked(0.0),Km(0.0),kdiff(0.0); //Diffusion coefficient [MJ/m2/K/d] + double md = 1.0; //Diffusion scaling factor in each layer + + //Water density in each layer as a function of water temperature. Reference: + // (Hostetler and Bartlein, 1990) double a0 = -1.954e-5; dens_e = DENSITY_WATER*(1.0+a0*pow(abs((T_old + ZERO_CELSIUS)-277.0),1.68)); dens_h = DENSITY_WATER*(1.0+a0*pow(abs((Ts_old + ZERO_CELSIUS)-277.0),1.68)); //Estimate layer thicknesses - double Ze = V_e_old/A_old; //Average epilimnion thickness [m] - double Zh = V_h_old/A_h_old; //Average hppolimnion thickness [m] - double Zc = (Ze+Zh)/2.0; //Characteristic length [m] + double Ze = V_e_old/A_old; //Average epilimnion thickness [m] + double Zh = V_h_old/A_h_old; //Average hypolimnion thickness [m] + double Zc = (Ze+Zh)/2.0; //Characteristic length [m] - //Molecular diffusion coefficient, Km can be calculated as per Community Land Model v5 (Lawrence et al, 2020) - Km = TC_WATER/HCP_WATER; //[MJ/m2/K/d] + //Molecular diffusion coefficient, Km can be calculated as per Community Land Model v5 + // (Lawrence et al, 2020) + Km = TC_WATER/HCP_WATER; // [m2/d] - //Estimate density gradient dp/dz over length Zc; ===> Brunt-Väisälä frequency across thermocline depth (Jennings et al. 2012) + //Estimate density gradient dp/dz over length Zc, and Brunt-Väisälä frequency across + // the thermocline depth (Jennings et al. 2012) if(dens_e FREEZING_TEMP){ - wstar = 0.0012*u2; //Surfaec friction velocity [m/s] - kstar = 6.6*pow(u2, -1.84) * sqrt(abs(sin(pHRU->GetLatRad()))); - Ri = 40.0*N2*pow(VON_KARMAN,2.0)*pow(Zc,2.0); - Ri /= pow(wstar,2.0); - Ri /= exp(-2*VON_KARMAN*Zc); - Ri += 1.0; - Ri = (-1 + sqrt(Ri))/20.0; - Ke = VON_KARMAN*wstar*Zc/(1.0 + 37.0*pow(Ri,2.0)) * exp(-VON_KARMAN*Zc); //[MJ/m2/K/s] - Ke *= SEC_PER_DAY; //[MJ/m2/K/d] + wstar = 0.0012*u2; //Surface friction velocity [m/s] + kstar = 6.6*pow(u2, -1.84) * sqrt(abs(sin(pHRU->GetLatRad()))); + Ri = 40.0*N2*pow(VON_KARMAN,2.0)*pow(zm,2.0); + Ri /= pow(wstar,2.0); + Ri /= exp(-2*kstar*zm); + Ri += 1.0; + Ri = (-1 + sqrt(Ri))/20.0; // Richardson number + if(Ri < 1e+20){ + Ke = VON_KARMAN*wstar*zm/(1.0 + 37.0*pow(Ri,2.0)) * exp(-kstar*zm); //[m2/s] + Ke *= SEC_PER_DAY; //[m2/d] + } } - //Enhanced diffusivity intended to represent unresolved mixing processes; ====> Community Land Model v5 (Lawrence et al, 2020) + //Enhanced diffusivity intended to represent unresolved mixing processes if (N2 >= 7.5e-5){ - Ked = 1.04e-8 * pow(N2,-0.43)*SEC_PER_DAY; //[MJ/m2/K/d] + Ked = 1.04e-8 * pow(N2,-0.43)*SEC_PER_DAY; //[m2/d] } //Increase the overall diffusivity for large lakes, intended to represent - //3-dimensional mixing processes such as cased by horizontal temperature; Gradient can be calculated as per Community Land Model v5 (Lawrence et al, 2020) + //3-dimensional mixing processes such as caused by horizontal temperature + // gradient. Calculated as per Community Land Model v5 (Lawrence et al, 2020) if ((Ze+Zh) >= 25.0){ - md = 10.0; //Diffusion scaling factor in each layer + md = 10.0; //Diffusion scaling factor } //Final diffusion coefficient - kdiff = md*(Ke+Ked+Km)*HCP_WATER; //[MJ/m2/K/d] - + //kdiff = md*(Ke+Ked+Km)*HCP_WATER; //[MJ/m2/K/d] + kdiff = threshMin(md*(Ke+Ked+Km)*HCP_WATER, 5.0, 0.0); //[MJ/m2/K/d] - //TEMPORARY OUPUT ====> To test which parameter is abnormal - //cout<<"Ze = "<0C water) E_guess +=dE; @@ -579,14 +629,12 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // iter++; //if (iter>2){cout<<"iter: "<GetHydroUnit(pRes->GetHRUIndex())->GetArea()*M2_PER_KM2; //double V_sed = pRes->GetLakebedThickness() * HRUarea; double V_new = pRes->GetStorage(); - double V_h_new = pRes->GetHypolimnionStorage(); + double V_h_new = pRes->GetHypolimnionStorage(); double V_e_new = V_new-V_h_new; mult = 1.0 / HRUarea; @@ -1324,18 +1379,21 @@ void CEnthalpyModel::WriteMinorOutput(const optStruct& Options,const time_struct double Qin=_pModel->GetSubBasin(p)->GetOutflowArray()[_pModel->GetSubBasin(p)->GetNumSegments()-1]*SEC_PER_DAY; //[m3/d] - GetEnergyLossesFromLake(p,Q_sens,Q_cond,Q_lat,Q_sw_in,Q_lw_in,Q_lw_out,Q_rain); + GetEnergyLossesFromLake(p,Q_sens,Q_conv,Q_lat,Q_sw_in,Q_lw_in,Q_lw_out,Q_rain,Q_adv,Ri,N2,kdiff); double lakeTemp =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_e_new); double sedTemp =ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_h_new); double inflowTemp=ConvertVolumetricEnthalpyToTemperature( Ein / Qin ); double pctFroz =ConvertVolumetricEnthalpyToIceContent (_aMres[p] / V_e_new); +// double stage2_mixing = _stage - _mixing_depth; _LAKEOUT << mult * Ein << "," << mult * Eout << ","; _LAKEOUT << mult * Q_rain << "," << mult * Q_sens << ","; - _LAKEOUT << mult * Q_cond << "," << mult * Q_lat << ","; + _LAKEOUT << mult * Q_conv << "," << mult * Q_adv << ","; + _LAKEOUT << mult * Ri << "," << N2 << "," << kdiff << ","; + _LAKEOUT << mult * Q_lat << ","; _LAKEOUT << mult * Q_sw_in << "," << mult * Q_lw_in << "," << mult * Q_lw_out << ","; - _LAKEOUT << mult * _aMres[p] << ","; + _LAKEOUT << mult * _aMres[p] << "," << mult * _aMsed[p] << ","; if (V_e_new!=0){_LAKEOUT << lakeTemp << ",";}else{_LAKEOUT<<",";} if (V_h_new!=0){_LAKEOUT << sedTemp << ",";}else{_LAKEOUT<<",";} if (Qin !=0){_LAKEOUT << inflowTemp << ",";}else{_LAKEOUT<<",";} diff --git a/src/EnergyTransport.h b/src/EnergyTransport.h index bd7af9f7..8cd074f7 100644 --- a/src/EnergyTransport.h +++ b/src/EnergyTransport.h @@ -60,7 +60,7 @@ class CEnthalpyModel :public CConstituentModel double GetEnergyLossesInTransit(const int p,double &Q_sens,double &Q_GW) const; double GetEnergyLossesFromReach(const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_GW,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out,double &Q_fric, double &Tave) const; - double GetEnergyLossesFromLake (const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out, double &Q_rain) const; + double GetEnergyLossesFromLake (const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out, double &Q_rain, double &Q_adv, double &Ri, double &N2, double &kdiff) const; double GetOutflowIceFraction (const int p) const; double GetAvgLatentHeatFlux () const; diff --git a/src/RavenInclude.h b/src/RavenInclude.h index 97d48c7c..fd52543c 100644 --- a/src/RavenInclude.h +++ b/src/RavenInclude.h @@ -179,7 +179,7 @@ const double TC_CLAY =0.216; const double TC_ORGANIC =0.0216; ///< [MJ/m/d/K] Thermal conductivity of organic matter (0.25 W/m/K) const double TC_AIR =0.00199; ///< [MJ/m/d/K] Thermal conductivity of air (0.023 W/m/K) -const double COM_WATER =4.58e-10; ///< [1/Pa] Compressiblity of Water +const double COM_WATER =4.58e-10; ///< [1/Pa] Compressiblity of Water const double COM_ICE =4.58e-10; ///< [1/Pa] Compressiblity of Ice const double HCP_WATER =4.186; ///< [MJ/m3/K] Volumetric Heat Capacity of Water const double HCP_ICE =1.938; ///< [MJ/m3/K] Volumetric Heat Capacity of Ice diff --git a/src/Reservoir.cpp b/src/Reservoir.cpp index 08bdf0eb..9e6dfcef 100755 --- a/src/Reservoir.cpp +++ b/src/Reservoir.cpp @@ -20,7 +20,7 @@ void CReservoir::BaseConstructor(const string Name,const long SubID) _pDownSB=NULL; _Q_dn_old = 0.0; _Q_up_old = 0.0; - _mixing_depth=5.0; + _mixing_depth=15.0; _lakebed_thick=1.0; _lakebed_cond =0.0; _lake_convcoeff=2.0; @@ -449,6 +449,11 @@ double CReservoir::GetSillElevation(const int nn) const } return _crest_ht+weir_adj; } +////////////////////////////////////////////////////////////////// +/// \returns mixing depth [m] +// +double CReservoir::GetMixingDepth () const { return _mixing_depth; } + ////////////////////////////////////////////////////////////////// /// \returns current surface area [m2] // diff --git a/src/Reservoir.h b/src/Reservoir.h index 8ae38567..5c231bab 100755 --- a/src/Reservoir.h +++ b/src/Reservoir.h @@ -206,6 +206,7 @@ class CReservoir double GetMinStage (const int nn) const;//[m] double GetMaxStage (const int nn) const;//[m] double GetSillElevation (const int nn) const;//[m] + double GetMixingDepth () const; //[m] int GetNumWaterDemands () const; CDemand *GetWaterDemandObj (const int ii) const; From 20b9596106472b0609bc330a3cc96e70e2bb2190 Mon Sep 17 00:00:00 2001 From: Markus Schnorbus Date: Thu, 28 Nov 2024 12:48:02 -0800 Subject: [PATCH 7/9] Mixing depth for each reservoir is now estimated as a*Fe^b, where the fetch, Fe, is estimated as the square root of maximum reservoir surface area (varies by reservoir) and the coefficient, a (mix_depth_coef), and exponent, b (mix_depth_expn), are provided as global parameters (MIX_DEPTH_COEF and MIX_DEPTH_EXPN, respectivley). --- src/EnergyTransport.cpp | 5 ++-- src/GlobalParams.cpp | 29 +++++++++++++++++++++ src/ParseHRUFile.cpp | 10 ++++---- src/ParsePropertyFile.cpp | 6 +++++ src/Properties.h | 3 +++ src/Reservoir.cpp | 54 +++++++++++++++++++++++++++++---------- src/Reservoir.h | 11 +++++--- 7 files changed, 92 insertions(+), 26 deletions(-) diff --git a/src/EnergyTransport.cpp b/src/EnergyTransport.cpp index 35a750bc..d3605ecf 100755 --- a/src/EnergyTransport.cpp +++ b/src/EnergyTransport.cpp @@ -375,7 +375,6 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub //kdiff = md*(Ke+Ked+Km)*HCP_WATER; //[MJ/m2/K/d] kdiff = threshMin(md*(Ke+Ked+Km)*HCP_WATER, 5.0, 0.0); //[MJ/m2/K/d] - Q_sens =hstar * A_avg * (temp_air -0.5*(T_new+T_old)); Q_conv =kdiff * A_h_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); Q_sw_in =(SW )*A_avg; @@ -459,7 +458,7 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // hstar =pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] // ksed =pRes->GetLakebedConductivity()/0.5/ pRes->GetLakebedThickness();//[MJ/m2/d/K] } - + ///////////// Option 1 to calculate kdiff //// water density in each layer as a function of water temperature (Chapra, 2008) //double dens_e,dens_h,vdiff,kdiff; @@ -1390,7 +1389,7 @@ void CEnthalpyModel::WriteMinorOutput(const optStruct& Options,const time_struct _LAKEOUT << mult * Ein << "," << mult * Eout << ","; _LAKEOUT << mult * Q_rain << "," << mult * Q_sens << ","; _LAKEOUT << mult * Q_conv << "," << mult * Q_adv << ","; - _LAKEOUT << mult * Ri << "," << N2 << "," << kdiff << ","; + _LAKEOUT << Ri << "," << N2 << "," << kdiff << ","; _LAKEOUT << mult * Q_lat << ","; _LAKEOUT << mult * Q_sw_in << "," << mult * Q_lw_in << "," << mult * Q_lw_out << ","; _LAKEOUT << mult * _aMres[p] << "," << mult * _aMsed[p] << ","; diff --git a/src/GlobalParams.cpp b/src/GlobalParams.cpp index 89c6ebfd..eb208999 100644 --- a/src/GlobalParams.cpp +++ b/src/GlobalParams.cpp @@ -156,6 +156,23 @@ void CGlobalParams::AutoCalculateGlobalParams(const global_struct &Gtmp, const g { G.GAMMA_SCALE_multiplier=1.0;//default=nothing (no warning needed) } + + autocalc=SetCalculableValue(G.mix_depth_coef,Gtmp.mix_depth_coef,Gtemplate.mix_depth_coef); + if (autocalc) + { + G.mix_depth_coef=1.1; + warn="The required global parameter MIX_DEPTH_COEF was autogenerated with value "+to_string(G.mix_depth_coef); + if (chatty){WriteAdvisory(warn,false);} + } + + autocalc=SetCalculableValue(G.mix_depth_expn,Gtmp.mix_depth_expn,Gtemplate.mix_depth_expn); + if (autocalc) + { + G.mix_depth_expn=5.3; + warn="The required global parameter MIX_DEPTH_EXPN was autogenerated with value "+to_string(G.mix_depth_expn); + if (chatty){WriteAdvisory(warn,false);} + } + autocalc=SetCalculableValue(G.reference_flow_mult,Gtmp.reference_flow_mult,Gtemplate.reference_flow_mult); if (autocalc) { @@ -417,6 +434,10 @@ void CGlobalParams::InitializeGlobalParameters(global_struct &g, bool is_templat g.TIME_TO_PEAK_multiplier =DefaultParameterValue(is_template,true); g.GAMMA_SHAPE_multiplier =DefaultParameterValue(is_template,true); g.GAMMA_SCALE_multiplier =DefaultParameterValue(is_template,true); + + g.mix_depth_coef =DefaultParameterValue(is_template,true); + g.mix_depth_expn =DefaultParameterValue(is_template,true); + //model-specific parameters g.avg_annual_snow =DefaultParameterValue(is_template,false); @@ -564,6 +585,10 @@ void CGlobalParams::SetGlobalProperty (global_struct &G, else if (!name.compare("HBVEC_LAPSE_RATE" )){G.HBVEC_lapse_rate=value; } else if (!name.compare("HBVEC_LAPSE_UPPER" )){G.HBVEC_lapse_upper=value; } else if (!name.compare("HBVEC_LAPSE_ELEV" )){G.HBVEC_lapse_elev=value; } + + else if (!name.compare("MIX_DEPTH_COEF" )){G.mix_depth_coef=value;} + else if (!name.compare("MIX_DEPTH_EXPN" )){G.mix_depth_expn=value;} + else{ WriteWarning("CGlobalParams::SetGlobalProperty: Unrecognized/invalid global parameter name ("+name+") in .rvp file",false); @@ -687,6 +712,10 @@ double CGlobalParams::GetGlobalProperty(const global_struct &G, string param_na else if (!name.compare("HBVEC_LAPSE_RATE" )){return G.HBVEC_lapse_rate; } else if (!name.compare("HBVEC_LAPSE_UPPER" )){return G.HBVEC_lapse_upper; } else if (!name.compare("HBVEC_LAPSE_ELEV" )){return G.HBVEC_lapse_elev; } + + else if (!name.compare("MIX_DEPTH_COEF" )){return G.mix_depth_coef;} + else if (!name.compare("MIX_DEPTH_EXPN" )){return G.mix_depth_expn;} + else{ if (strict){ string msg="CGlobalParams::GetParameter: Unrecognized/invalid global parameter name in .rvp file: "+name; diff --git a/src/ParseHRUFile.cpp b/src/ParseHRUFile.cpp index 7c12d1a2..12bdd1ef 100644 --- a/src/ParseHRUFile.cpp +++ b/src/ParseHRUFile.cpp @@ -412,7 +412,7 @@ bool ParseHRUPropsFile(CModel *&pModel, const optStruct &Options, bool terrain_r //in*=CGlobalParams::GetParameter("TOC_MULTIPLIER"); // use to be this; but has its own multiplier now; maybe set default to TOC_MULTIPLIER?? in *= pModel->GetGlobalParams()->GetParameter("TIME_TO_PEAK_MULTIPLIER"); } - if(!aParamStrings[i].compare("GAMMA_SHAPE") && (in!=AUTO_COMPUTE) && (in!=USE_TEMPLATE_VALUE)){ + if(!aParamStrings[i].compare("GAMMA_SHAPE") && (in!=AUTO_COMPUTE) && (in!=USE_TEMPLATE_VALUE)){ in *= pModel->GetGlobalParams()->GetParameter("GAMMA_SHAPE_MULTIPLIER"); } if(!aParamStrings[i].compare("GAMMA_SCALE") && (in!=AUTO_COMPUTE) && (in!=USE_TEMPLATE_VALUE)){ @@ -1882,15 +1882,15 @@ CReservoir *ReservoirParse(CParser *p,string name,const CModel *pModel,long long //------------------------------------------------------------------------------------ if((type==CURVE_POWERLAW) || (type==CURVE_LINEAR)) { - pRes=new CReservoir(name,SBID,a_V,b_V,a_Q,b_Q,a_A,b_A,crestht,max_depth); + pRes=new CReservoir(name,SBID,a_V,b_V,a_Q,b_Q,a_A,b_A,crestht,max_depth,pModel); } else if(type==CURVE_DATA) { - pRes=new CReservoir(name,SBID,aQ_ht,aQ,aQund,aA,aV,NQ);//presumes aQ_ht=aV_ht=aA_ht; NA=NV=NQ + pRes=new CReservoir(name,SBID,aQ_ht,aQ,aQund,aA,aV,NQ,pModel);//presumes aQ_ht=aV_ht=aA_ht; NA=NV=NQ } else if(type==CURVE_VARYING) { - pRes=new CReservoir(name,SBID,nDates,aDates,aQ_ht,aQQ,aQund,aA,aV,NQ);//presumes aQ_ht=aV_ht=aA_ht; NA=NV=NQ + pRes=new CReservoir(name,SBID,nDates,aDates,aQ_ht,aQQ,aQund,aA,aV,NQ,pModel);//presumes aQ_ht=aV_ht=aA_ht; NA=NV=NQ } else if(type==CURVE_LAKE) { @@ -1907,7 +1907,7 @@ CReservoir *ReservoirParse(CParser *p,string name,const CModel *pModel,long long } } - pRes=new CReservoir(name,SBID,weircoeff,cwidth,crestht,lakearea,max_depth); + pRes=new CReservoir(name,SBID,weircoeff,cwidth,crestht,lakearea,max_depth,pModel); } else { diff --git a/src/ParsePropertyFile.cpp b/src/ParsePropertyFile.cpp index 1ee704cd..13604f3f 100644 --- a/src/ParsePropertyFile.cpp +++ b/src/ParsePropertyFile.cpp @@ -167,6 +167,12 @@ bool ParseClassPropertiesFile(CModel *&pModel, aP [0]="GAMMA_SCALE_MULTIPLIER"; aPC[0]=CLASS_GLOBAL; AddToMasterParamList(aPmaster, aPCmaster,nPmaster, aP, aPC, 1); + + aP [0]="MIX_DEPTH_COEF"; aPC[0]=CLASS_GLOBAL; + AddToMasterParamList(aPmaster, aPCmaster,nPmaster, aP, aPC, 1); + + aP [0]="MIX_DEPTH_EXPN"; aPC[0]=CLASS_GLOBAL; + AddToMasterParamList(aPmaster, aPCmaster,nPmaster, aP, aPC, 1); //Throw warning if NULL Terrain but terrain parameter is needed terrain_required=false; diff --git a/src/Properties.h b/src/Properties.h index 52a210d2..71b0f7fb 100644 --- a/src/Properties.h +++ b/src/Properties.h @@ -396,6 +396,7 @@ struct global_struct double avg_annual_snow; ///< [mm] avg annual snow as SWE double avg_annual_runoff; ///< [mm] avg annual runoff from basin + double reference_flow_mult; ///< [-] multiplier for reference flow relative to annual mean flow ///< (10 for flood modelling, 0.5 for low-flow conditions) @@ -426,6 +427,8 @@ struct global_struct double TIME_TO_PEAK_multiplier; ///< [0..1+] time to peak multiplier double GAMMA_SHAPE_multiplier; ///< [0..1+] Gamma shape multiplier double GAMMA_SCALE_multiplier; ///< [0..1+] Gamma scale multiplier + double mix_depth_coef; ///< [m] Reservoir mixing depth coefficient + double mix_depth_expn; ///< [-] Reservoir mixing depth exponent double assimilation_fact; ///< [0..1] assimilation factor (0=no assimilation to 1= full replacement) double assim_upstream_decay; ///< [1/km] assimilation upstream decay factor (0= overrides everything upstream, large- observation influence decays quickly with distance from gauge) [~0.01] double assim_time_decay; ///< [1/d] assimilation temporal decay factor (0=diminishes in future, 0.1 - diminshes in 3 days) [0.2] diff --git a/src/Reservoir.cpp b/src/Reservoir.cpp index 9e6dfcef..010e2585 100755 --- a/src/Reservoir.cpp +++ b/src/Reservoir.cpp @@ -11,7 +11,7 @@ /// \param SubID [in] subbasin ID /// \param typ [in] reservoir type /// \details needed because versions of c++ prior to v11 don't necessaarily support delegating constructors -// +// void CReservoir::BaseConstructor(const string Name,const long SubID) { _name=Name; @@ -20,7 +20,7 @@ void CReservoir::BaseConstructor(const string Name,const long SubID) _pDownSB=NULL; _Q_dn_old = 0.0; _Q_up_old = 0.0; - _mixing_depth=15.0; + _mixing_depth=0.0; _lakebed_thick=1.0; _lakebed_cond =0.0; _lake_convcoeff=2.0; @@ -32,7 +32,7 @@ void CReservoir::BaseConstructor(const string Name,const long SubID) _Qout =0.0; _Qout_last =0.0; _MB_losses =0.0; - _AET =0.0; + _AET =0.0; _Precip =0.0; _GW_seepage=0.0; _aQstruct=NULL; @@ -87,6 +87,7 @@ void CReservoir::BaseConstructor(const string Name,const long SubID) _nControlStructures=0; _pControlStructures=NULL; + _seepage_const=0; _local_GW_head=0.0; @@ -119,15 +120,17 @@ CReservoir::CReservoir(const string Name,const long SubID) /// \param b_V [in] power law exponent for volume rating curve [-] /// \param a_Q [in] power law coefficient for discharge rating curve [m3/s*m^-b_Q] /// \param b_Q [in] power law exponent for discharge rating curve [-] -/// \param a_A[in] surface area of reservoir/lake when stage at absolute crest height [m2] +/// \param a_A [in] surface area of reservoir/lake when stage at absolute crest height [m2] /// \param b_A [in] power law exponent for area rating curve [-] (0 for prismatic reservoir) // + CReservoir::CReservoir(const string Name, const long SBID, const double a_V, const double b_V, const double a_Q, const double b_Q, const double a_A, const double b_A, - const double crestht,const double depth) + const double crestht,const double depth, + const CModelABC* pModel) { BaseConstructor(Name,SBID); @@ -141,6 +144,10 @@ CReservoir::CReservoir(const string Name, const long SBID, double aA=a_A; double bA=b_A; + + double COEFF_RER_A = pModel->GetGlobalParams()->GetParams()->mix_depth_coef; + double COEFF_RER_B = pModel->GetGlobalParams()->GetParams()->mix_depth_expn; + if(aA==AUTO_COMPUTE) {aA=a_V/depth*b_V;} if(bA==AUTO_COMPUTE) {bA=b_V-1.0; } @@ -149,6 +156,7 @@ CReservoir::CReservoir(const string Name, const long SBID, _min_stage =_crest_ht-depth; _max_stage =_crest_ht+6.0; // a postive value relative to _crest_ht dh=(_max_stage-_min_stage)/(double)(_Np-1); + _mixing_depth = COEFF_RER_A*pow(sqrt(aA),COEFF_RER_B); for(int i=0;i<_Np;i++) { @@ -169,13 +177,13 @@ CReservoir::CReservoir(const string Name, const long SBID, /// \param SubID [in] subbasin ID /// \param a_ht[] [in] array of reservoir stages [size: nPoints] /// \param a_Q[] [in] array of reservoir discharges [size: nPoints] -/// \param a_A[] [in] array of reservoir volumes [size: nPoints] -/// \param a_V[] [in] array of reservoir areas [size: nPoints] +/// \param a_A[] [in] array of reservoir areas [size: nPoints] +/// \param a_V[] [in] array of reservoir volumes [size: nPoints] // CReservoir::CReservoir(const string Name, const long SubID, const double *a_ht, const double *a_Q, const double *a_Qund,const double *a_A, const double *a_V, - const int nPoints) + const int nPoints,const CModelABC* pModel) { BaseConstructor(Name,SubID); @@ -191,6 +199,8 @@ CReservoir::CReservoir(const string Name, const long SubID, _aArea =new double [_Np]; _aVolume=new double [_Np]; ExitGracefullyIf(_aVolume==NULL,"CReservoir::constructor (2)",OUT_OF_MEMORY); + double COEFF_RER_A = pModel->GetGlobalParams()->GetParams()->mix_depth_coef; + double COEFF_RER_B = pModel->GetGlobalParams()->GetParams()->mix_depth_expn; string warn; for (int i=0;i<_Np;i++) @@ -225,6 +235,7 @@ CReservoir::CReservoir(const string Name, const long SubID, } } _max_capacity=_aVolume[_Np-1]; + _mixing_depth = COEFF_RER_A*pow(sqrt(_aArea[_Np-1]),COEFF_RER_B); } ////////////////////////////////////////////////////////////////// @@ -245,7 +256,8 @@ CReservoir::CReservoir(const string Name, const long SubID, const double *a_Qund, const double *a_A, const double *a_V, - const int nPoints) + const int nPoints, + const CModelABC* pModel) { BaseConstructor(Name,SubID); @@ -270,6 +282,9 @@ CReservoir::CReservoir(const string Name, const long SubID, _aVolume=new double [_Np]; ExitGracefullyIf(_aVolume==NULL,"CReservoir::constructor (2)",OUT_OF_MEMORY); + double COEFF_RER_A = pModel->GetGlobalParams()->GetParams()->mix_depth_coef; + double COEFF_RER_B = pModel->GetGlobalParams()->GetParams()->mix_depth_expn; + string warn; for (int i=0;i<_Np;i++) { @@ -302,6 +317,7 @@ CReservoir::CReservoir(const string Name, const long SubID, } } _max_capacity=_aVolume[_Np-1]; + _mixing_depth = COEFF_RER_A*pow(sqrt(_aArea[_Np-1]),COEFF_RER_B); } ////////////////////////////////////////////////////////////////// /// \brief Constructor for prismatic lake reservoir controlled by weir coefficient @@ -318,7 +334,8 @@ CReservoir::CReservoir(const string Name, const double crestw, const double crestht, const double A, - const double depth) + const double depth, + const CModelABC* pModel) { BaseConstructor(Name,SubID); @@ -340,7 +357,10 @@ CReservoir::CReservoir(const string Name, _aArea =new double [_Np]; _aVolume=new double [_Np]; ExitGracefullyIf(_aVolume==NULL,"CReservoir::constructor (4)",OUT_OF_MEMORY); - + + double COEFF_RER_A = pModel->GetGlobalParams()->GetParams()->mix_depth_coef; + double COEFF_RER_B = pModel->GetGlobalParams()->GetParams()->mix_depth_expn; + double dh; string warn; dh=(_max_stage-_crest_ht)/(_Np-2); //spacing = 0.05m @@ -348,6 +368,7 @@ CReservoir::CReservoir(const string Name, _aQ [0]=0.0; _aQunder[0]=0.0; _aArea [0]=A; + _aVolume[0]=0.0; for (int i=1;i<_Np;i++) // - Edited by KL to reference crest height properly { @@ -358,6 +379,7 @@ CReservoir::CReservoir(const string Name, _aVolume[i]=A*(_aStage[i]-_min_stage); } _max_capacity=_aVolume[_Np-1]; + _mixing_depth = COEFF_RER_A*pow(sqrt(A),COEFF_RER_B); } ////////////////////////////////////////////////////////////////// @@ -400,6 +422,10 @@ CReservoir::~CReservoir() // long CReservoir::GetSubbasinID () const { return _SBID; } +/// \returns mixing depth [m] +// +double CReservoir::GetMixingDepth () const { return _mixing_depth;} + ////////////////////////////////////////////////////////////////// /// \returns reservoir storage [m3] // @@ -450,9 +476,6 @@ double CReservoir::GetSillElevation(const int nn) const return _crest_ht+weir_adj; } ////////////////////////////////////////////////////////////////// -/// \returns mixing depth [m] -// -double CReservoir::GetMixingDepth () const { return _mixing_depth; } ////////////////////////////////////////////////////////////////// /// \returns current surface area [m2] @@ -574,6 +597,7 @@ double CReservoir::GetCrestWidth() const { return _crest_width; } + ////////////////////////////////////////////////////////////////// /// \brief gets max capacity /// \returns capacity, in m3 @@ -958,6 +982,8 @@ void CReservoir::SetCrestWidth(const double& width) _crest_width=width; } } + + ////////////////////////////////////////////////////////////////// /// \brief sets max capacity /// \param capacity, in m3 diff --git a/src/Reservoir.h b/src/Reservoir.h index 5c231bab..d3d8f392 100755 --- a/src/Reservoir.h +++ b/src/Reservoir.h @@ -163,17 +163,20 @@ class CReservoir const double a_v, const double b_v, const double a_Q, const double b_Q, const double a_A, const double b_A, - const double crestht, const double depth); + const double crestht, const double depth, + const CModelABC* pModel); CReservoir(const string name, const long SubID, const double *a_ht, const double *a_Q, const double *aQ_und,const double *a_A, const double *a_V, - const int nPoints); + const int nPoints, + const CModelABC* pModel); CReservoir(const string Name, const long SubID, const int nDates, const int *aDates, const double *a_ht, double **a_QQ, const double *aQ_und, const double *a_A, const double *a_V, - const int nPoints); + const int nPoints, + const CModelABC* pModel); CReservoir(const string Name, const long SubID, const double weircoeff, //Lake constructor - const double crestw, const double crestht, const double A, const double depth); + const double crestw, const double crestht, const double A, const double depth,const CModelABC* pModel); ~CReservoir(); //Accessors From a0b306a9916308b6236023c4b79f958c4459286e Mon Sep 17 00:00:00 2001 From: Markus Schnorbus Date: Fri, 31 Jan 2025 11:39:08 -0800 Subject: [PATCH 8/9] Modifications to reservoir energy routing to more efficiently allow two options to estimate kdiff, added tracking of previous vertical flow rate, corrected solution for use of proper previous energy storage in epilimnion and hypolimnion, added estimate of dE and dEs in case numerical algoriyhm fails to converge, modified calaculation of lake energy balances for output. --- src/EnergyTransport.cpp | 364 +++++++++++++++++++++------------------- src/EnergyTransport.h | 2 +- src/Reservoir.cpp | 41 +++++ src/Reservoir.h | 11 +- 4 files changed, 239 insertions(+), 179 deletions(-) diff --git a/src/EnergyTransport.cpp b/src/EnergyTransport.cpp index 339348fd..9703b0b8 100755 --- a/src/EnergyTransport.cpp +++ b/src/EnergyTransport.cpp @@ -235,6 +235,7 @@ void CEnthalpyModel::SetBedTemperature(const int p,const double &val) { _aBedTemp[p]=val; } + ////////////////////////////////////////////////////////////////// /// \brief Set hyporheic layer index /// \param m [in] soil layer @@ -246,18 +247,31 @@ void CEnthalpyModel::SetHyporheicLayer(const int m) } _mHyporheic=m; } + ////////////////////////////////////////////////////////////////// /// \brief Calculates individual energy gain terms from lake/reservoir over current time step -/// \param p subbasin index -/// \param Q_sens [out] energy gain from convective/sensible heating [MJ/d] -/// \param Q_cond [out] energy gain from conductive exchange with bed [MJ/d] -/// \param Q_lat [out] energy gain from latent heat exchange [MJ/d] -/// \param Q_rad [out] energy gain from SW/LW radiation at surface [MJ/d] -/// \param Q_lw_out [out] (negative) energy gain from -/// \param Q_rain [out] energy gain from precip inputs [MJ/d] -/// \returns total energy lost from reach over current time step [MJ] +/// \param p [in] subbasin index +/// \param Q_sens [out] energy gain from convective/sensible heating [MJ/d] +/// \param Q_conv [out] energy gain from convective exchange between layers [MJ/d] +/// \param Q_lat [out] energy gain from latent heat exchange [MJ/d] +/// \param Q_sw_in [out] energy gain from incoming SW radiation at surface [MJ/d] +/// \param Q_lw_in [out] energy gain from incoming LW radiation at surface [MJ/d] +/// \param Q_lw_out[out] (negative) energy gain from LW radiation at surface [MJ/d] +/// \param Q_rain [out] energy gain from precip inputs [MJ/d] +/// \param Q_adv [out] energy gain advective heat exchange between layers [MJ/d] +/// \param kdiff [out] coefficient controlling convection between layers [MJ/m2/d/K] +/// \returns total energy lost from lake/reservoir over current time step [MJ] // -double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, double &Q_conv, double &Q_lat, double &Q_sw_in, double &Q_lw_in, double &Q_lw_out, double &Q_rain, double &Q_adv, double &Ri, double &N2, double &kdiff) const +double CEnthalpyModel::GetEnergyLossesFromLake(const int p, + double &Q_sens, + double &Q_conv, + double &Q_lat, + double &Q_sw_in, + double &Q_lw_in, + double &Q_lw_out, + double &Q_rain, + double &Q_adv, + double &kdiff) const { double tstep = _pModel->GetOptStruct()->timestep; @@ -267,22 +281,14 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); double V_new = pRes->GetStorage(); - double V_old = pRes->GetOldStorage(); - double Q_new = pRes->GetOutflowRate() * SEC_PER_DAY; - double Q_old = pRes->GetOldOutflowRate() * SEC_PER_DAY; double A_new = pRes->GetSurfaceArea(); - double A_old = pRes->GetOldSurfaceArea(); - double A_avg = 0.5 * (A_new + A_old); double V_h_new = pRes->GetHypolimnionStorage(); double V_h_old = pRes->GetOldHypolimnionStorage(); double V_e_new = V_new-V_h_new; - double V_e_old = V_old-V_h_old; double A_h_new = pRes->GetMixingArea(); - double A_h_old = pRes->GetOldMixingArea(); - double A_h_avg = 0.5 * (A_h_new + A_h_old); double zm = pRes->GetMixingDepth(); - double Q_dn_new(0),Q_dn_old(0),Q_up_new(0),Q_up_old(0); + double Q_dn_new(0),Q_up_new(0); double Q_vert = (V_h_new - V_h_old)/tstep; if (Q_vert> 0){ Q_dn_new = Q_vert; Q_up_new = 0.0; @@ -293,98 +299,30 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub double Acorr=1.0; double SW(0), LW(0), LW_in(0), temp_air(0), AET(0); double hstar(0), ksed(0), Vsed=0.001; - double u2(0); - double T_new =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_e_new); - double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p] / V_e_old); double Ts_new=ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_h_new ); - double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p] / V_h_old ); - if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input - Acorr =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area + if(pHRU!=NULL) { //otherwise only simulate advective mixing+rain input + Acorr =pHRU->GetArea()*M2_PER_KM2/A_new; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area temp_air =pHRU->GetForcingFunctions()->temp_ave; //[C] SW =pHRU->GetForcingFunctions()->SW_radia_net; //[MJ/m2/d] - not using canopy correction! LW_in =pHRU->GetForcingFunctions()->LW_incoming; //[MJ/m2/d] - LW =-STEFAN_BOLTZ*EMISS_WATER*0.5*(pow(T_new+ZERO_CELSIUS,4)+pow(T_old+ZERO_CELSIUS,4)); + LW =-STEFAN_BOLTZ*EMISS_WATER*pow(T_new+ZERO_CELSIUS,4); //[MJ/m2/d] AET =pRes->GetAET()*Acorr/ MM_PER_METER ; //[m/d] //*pHRU->GetArea()/A_avg - hstar =pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] - Vsed = pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; - ksed = pRes->GetLakebedConductivity() / 0.5 / pRes->GetLakebedThickness(); //[MJ/m2/d/K] + hstar =pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K] + Vsed =pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2; + ksed =pRes->GetLakebedConductivity() / 0.5 / pRes->GetLakebedThickness(); //[MJ/m2/d/K] } + kdiff = pRes->GetLakeLyrConvectionCoeff(); //[MJ/m2/d/K] - ///////////// Option 2 to calculate kdiff===> kdiff = md*(Ke+Ked+Km); - u2 = pHRU->GetForcingFunctions()->wind_vel; - - double kstar = 0.0; - double wstar = 0.0; //Surface friction velocity [m/s] - Ri = 0.0; //Richardson number - N2 = 5.0E-12; //Brunt-Väisälä frequency when dens_e >= dens_h [1/s2] - double dens_e, dens_h; //Water density in each layer [kg/m3] - double Ke(0.0),Ked(0.0),Km(0.0); //kdiff(0.0); //Diffusion coefficient [MJ/m2/K/d] - kdiff = 0.0; - double md = 1.0; //Diffusion scaling factor in each layer - - //Water density in each layer as a function of water temperature. Reference: - // (Hostetler and Bartlein, 1990) - double a0 = -1.954e-5; - dens_e = DENSITY_WATER*(1.0+a0*pow(abs((T_old + ZERO_CELSIUS)-277.0),1.68)); - dens_h = DENSITY_WATER*(1.0+a0*pow(abs((Ts_old + ZERO_CELSIUS)-277.0),1.68)); - - //Estimate layer thicknesses - double Ze = V_e_old/A_old; //Average epilimnion thickness [m] - double Zh = V_h_old/A_h_old; //Average hppolimnion thickness [m] - double Zc = (Ze+Zh)/2.0; //Characteristic length [m] - - //Molecular diffusion coefficient, Km can be calculated as per Community Land Model v5 - //(Lawrence et al, 2020) - Km = TC_WATER/HCP_WATER; //[MJ/m2/K/d] - - //Estimate density gradient dp/dz over length Zc; ===> Brunt-Väisälä frequency across - //the thermocline depth (Jennings et al. 2012) - if(dens_e FREEZING_TEMP){ - wstar = 0.0012*u2; //Surface friction velocity [m/s] - kstar = 6.6*pow(u2, -1.84) * sqrt(abs(sin(pHRU->GetLatRad()))); - Ri = 40.0*N2*pow(VON_KARMAN,2.0)*pow(zm,2.0); - Ri /= pow(wstar,2.0); - Ri /= exp(-2*kstar*zm); - Ri += 1.0; - Ri = (-1 + sqrt(Ri))/20.0; // Richardson number - if(Ri < 1e+20){ - Ke = VON_KARMAN*wstar*zm/(1.0 + 37.0*pow(Ri,2.0)) * exp(-kstar*zm); //[MJ/m2/K/s] - Ke *= SEC_PER_DAY; //[MJ/m2/K/d] - } - } - - //Enhanced diffusivity intended to represent unresolved mixing processes - if (N2 >= 7.5e-5){ - Ked = 1.04e-8 * pow(N2,-0.43)*SEC_PER_DAY; //[MJ/m2/K/d] - } - - //Increase the overall diffusivity for large lakes, intended to represent - //3-dimensional mixing processes such as cased by horizontal temperature; - // Gradient can be calculated as per Community Land Model v5 (Lawrence et al, 2020) - if ((Ze+Zh) >= 25.0){ - md = 10.0; //Diffusion scaling factor in each layer - } - - //Final diffusion coefficient - //kdiff = md*(Ke+Ked+Km)*HCP_WATER; //[MJ/m2/K/d] - kdiff = threshMin(md*(Ke+Ked+Km)*HCP_WATER, 5.0, 0.0); //[MJ/m2/K/d] - - Q_sens =hstar * A_avg * (temp_air -0.5*(T_new+T_old)); - Q_conv =kdiff * A_h_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old)); - Q_sw_in =(SW )*A_avg; - Q_lw_in =(LW_in )*A_avg; - Q_lw_out=(LW )*A_avg; - Q_lat =-(AET * DENSITY_WATER * LH_VAPOR)*A_avg; - Q_rain =_aMresRain[p]; - Q_adv = (Q_up_new*_aMsed[p]/V_e_new - Q_dn_new*_aMres[p]/V_h_new) * tstep; + Q_sens =hstar * A_new * (temp_air - T_new); //[MJ/d] + Q_conv =kdiff * A_h_new * (Ts_new - T_new); //[MJ/d] + Q_sw_in =(SW )*A_new; //[MJ/d] + Q_lw_in =(LW_in )*A_new; //[MJ/d] + Q_lw_out=(LW )*A_new; //[MJ/d] + Q_lat =-(AET * DENSITY_WATER * LH_VAPOR)*A_new; //[MJ/d] + Q_rain =_aMresRain[p]; //[MJ/d] + Q_adv = (Q_up_new*_aMsed[p]/V_e_new - Q_dn_new*_aMres[p]/V_h_new); //[MJ/d] //return -(Q_sens + Q_cond + Q_sw_in + Q_lw_in + Q_lw_out + Q_lat + Q_rain) * tstep; //[MJ] return -(Q_sens + Q_sw_in + Q_lw_in + Q_lw_out + Q_lat + Q_rain) * tstep; //[MJ] @@ -422,6 +360,8 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // double V_e_old = V_old-V_h_old; double Q_new = pRes->GetOutflowRate() * SEC_PER_DAY; double Q_old = pRes->GetOldOutflowRate() * SEC_PER_DAY; + double Q_dn_old = pRes->GetOldDownwardFlowRate(); + double Q_up_old = pRes->GetOldUpwardFlowRate(); double A_new = pRes->GetSurfaceArea(); double A_old = pRes->GetOldSurfaceArea(); double A_avg = 0.5 * (A_new + A_old); @@ -432,7 +372,7 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // double tstep=Options.timestep; - double Q_dn_new(0),Q_dn_old(0),Q_up_new(0),Q_up_old(0); + double Q_dn_new(0), Q_up_new(0); double Q_vert = (V_h_new - V_h_old)/tstep; if (Q_vert> 0){ Q_dn_new = Q_vert; Q_up_new = 0.0; @@ -442,12 +382,13 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // CHydroUnit* pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex()); double Acorr=1.0; - double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p]/V_e_old); - double Ts_old =ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p]/V_h_old); + double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres[p]/V_e_old); + double Ts_old =ConvertVolumetricEnthalpyToTemperature(_aMsed[p]/V_h_old); double SW(0), LW(0), LW_in(0), temp_air(0), AET(0); double hstar(0); // ksed=0.0, Vsed=0.001; - double u2(0); + double dens_e(0), dens_h(0); //Water density in each layer [kg/m3] + double N2 = 5.0E-12; //Brunt-Väisälä frequency when dens_e >= dens_h [1/s2] if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input Acorr =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area @@ -461,62 +402,44 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // // ksed =pRes->GetLakebedConductivity()/0.5/ pRes->GetLakebedThickness();//[MJ/m2/d/K] } -///////////// Option 1 to calculate kdiff -//// water density in each layer as a function of water temperature (Chapra, 2008) -//double dens_e,dens_h,vdiff,kdiff; -//double a0= 999.842594; -//double a1= 6.793952e-2; -//double a2=-9.095290e-3; -//double a3= 1.001685e-4; -//double a4=-1.120083e-6; -//double a5= 6.536332e-9; -//dens_e = a0+a1*T_old + a2*pow(T_old,2) + a3*pow(T_old,3) + a4*pow(T_old,4) + a5*pow(T_old,5); -//dens_h = a0+a1*Ts_old + a2*pow(Ts_old,2) + a3*pow(Ts_old,3) + a4*pow(Ts_old,4) + a5*pow(Ts_old,5); -//// Brunt-Väisälä frequency at thermocline depth (assume dz = 1.0m) (Jennings et al. 2012) -//double N2=1.0e-12; // value when dens_e >= dens_h -//if(dens_e kdiff = md*(Ke+Ked+Km); - u2 = pHRU->GetForcingFunctions()->wind_vel; - - double kstar = 0.0; - double wstar = 0.0; //Surface friction velocity [m/s] - double Ri = 0.0; //Richardson number - double N2 = 5.0E-12; //Brunt-Väisälä frequency when dens_e >= dens_h [1/s2] - double dens_e, dens_h; //Water density in each layer [kg/m3] - double Ke(0.0),Ked(0.0),Km(0.0),kdiff(0.0); //Diffusion coefficient [MJ/m2/K/d] - double md = 1.0; //Diffusion scaling factor in each layer - //Water density in each layer as a function of water temperature. Reference: // (Hostetler and Bartlein, 1990) double a0 = -1.954e-5; dens_e = DENSITY_WATER*(1.0+a0*pow(abs((T_old + ZERO_CELSIUS)-277.0),1.68)); dens_h = DENSITY_WATER*(1.0+a0*pow(abs((Ts_old + ZERO_CELSIUS)-277.0),1.68)); - //Estimate layer thicknesses + //Estimate layer thicknesses double Ze = V_e_old/A_old; //Average epilimnion thickness [m] double Zh = V_h_old/A_h_old; //Average hypolimnion thickness [m] double Zc = (Ze+Zh)/2.0; //Characteristic length [m] - //Molecular diffusion coefficient, Km can be calculated as per Community Land Model v5 - // (Lawrence et al, 2020) - Km = TC_WATER/HCP_WATER; // [m2/d] - //Estimate density gradient dp/dz over length Zc, and Brunt-Väisälä frequency across // the thermocline depth (Jennings et al. 2012) if(dens_e Vertical convection coefficient as function of N2 (Quay et el. 1980; Figure 10) + //double vdiff(0.0), kdiff(0.0); + //double m = -0.65; + //double b = -3.1; + //vdiff = pow(10,m*log10(N2)+b)/100/100; + //kdiff = vdiff*HCP_WATER*SEC_PER_DAY; + ///////////////// + + /// Option 2 --> kdiff = md*(Ke+Ked+Km) as per Community Land Model v5 (Lawrence et al, 2020) + double u2 = pHRU->GetForcingFunctions()->wind_vel; + double kstar = 0.0; + double wstar = 0.0; //Surface friction velocity [m/s] + double Ri = 0.0; //Richardson number + double Ke(0.0),Ked(0.0),Km(0.0),kdiff(0.0); //Diffusion coefficient [MJ/m2/K/d] + double md = 1.0; //Diffusion scaling factor in each layer + + //Molecular diffusion coefficient, Km can be calculated as per Community Land Model v5 + // (Lawrence et al, 2020) + Km = TC_WATER/HCP_WATER; // [m2/d] + //Wind-driven eddy diffusion coefficient, Ke, calculated as per per Community Land //Model v5 (Lawrence et al, 2020) if(T_old > FREEZING_TEMP){ @@ -542,12 +465,16 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // //3-dimensional mixing processes such as caused by horizontal temperature // gradient. Calculated as per Community Land Model v5 (Lawrence et al, 2020) if ((Ze+Zh) >= 25.0){ - md = 10.0; //Diffusion scaling factor + //md = 10.0; //Diffusion scaling factor + md = 1.0; } //Final diffusion coefficient - //kdiff = md*(Ke+Ked+Km)*HCP_WATER; //[MJ/m2/K/d] kdiff = threshMin(md*(Ke+Ked+Km)*HCP_WATER, 5.0, 0.0); //[MJ/m2/K/d] + //kdiff = md*(Ke+Ked+Km)*HCP_WATER; //[MJ/m2/K/d] + //////////////////////////////// + + pRes->SetLakeLyrConvectionCoeff(kdiff); // N-R solution of Crank-nicolson problem as set of two non-linear algebraic equations [A][E]=[B] // ----------------------------------------------------------------------------------------- @@ -567,32 +494,32 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // } //Constants matrix - B[0] = 0.5*(aMout_new[nSegments-1]+_aMout[p][nSegments-1])*tstep; // inflow [MJ] - B[0]+=(1.0-0.5*tstep/V_e_old*Q_old)*_aMres_last[p]; // outflow [MJ] - B[0]+=_aMresRain[p]*tstep; // rainfall inputs [MJ] - B[0]+= A_avg*(SW+LW_in)*tstep; // net incoming radiation [MJ] - B[0]+= A_avg*(LW )*tstep; // net outgoing radiation [MJ] - B[0]-= A_avg*(AET*DENSITY_WATER*LH_VAPOR)*tstep; // latent heat [MJ] - B[0]+= A_avg*hstar*(temp_air-0.5*T_old)*tstep; // sensible heat exhange [MJ] - B[0]+=0.5*kdiff*A_h_avg*(Ts_old-T_old)*tstep; // diffusive heat exchange [MJ] - B[0]-=0.5*tstep/V_e_old*Q_dn_old*_aMres_last[p]*tstep; // downward advection [MJ] - B[0]+=0.5*Q_up_old/V_h_old*_aMsed_last[p]*tstep; // upward advection [MJ] + B[0] = 0.5*(aMout_new[nSegments-1]+_aMout[p][nSegments-1])*tstep; // inflow [MJ] + B[0]+=(1.0-0.5*tstep/V_e_old*Q_old)*_aMres[p]; // outflow [MJ] + B[0]+=_aMresRain[p]*tstep; // rainfall inputs [MJ] + B[0]+= A_avg*(SW+LW_in)*tstep; // net incoming radiation [MJ] + B[0]+= A_avg*(LW )*tstep; // net outgoing radiation [MJ] + B[0]-= A_avg*(AET*DENSITY_WATER*LH_VAPOR)*tstep; // latent heat [MJ] + B[0]+= A_avg*hstar*(temp_air-0.5*T_old)*tstep; // sensible heat exhange [MJ] + B[0]+=0.5*kdiff*A_h_avg*(Ts_old-T_old)*tstep; // diffusive heat exchange [MJ] + B[0]-=0.5*tstep/V_e_old*Q_dn_old*_aMres[p]*tstep; // downward advection [MJ] + B[0]+=0.5*Q_up_old/V_h_old*_aMsed[p]*tstep; // upward advection [MJ] - B[1] =0.5*kdiff*A_h_avg*(T_old-Ts_old)*tstep; // diffusive heat exchange [MJ] - B[1]+=0.5/V_e_old*Q_dn_old*_aMres_last[p]*tstep; // downward advection [MJ] - B[1]+=(1-0.5*tstep/V_h_old*Q_up_old)*_aMsed_last[p]; // upward advection [MJ] + B[1] =0.5*kdiff*A_h_avg*(T_old-Ts_old)*tstep; // diffusive heat exchange [MJ] + B[1]+=0.5/V_e_old*Q_dn_old*_aMres[p]*tstep; // downward advection [MJ] + B[1]+=(1-0.5*tstep/V_h_old*Q_up_old)*_aMsed[p]; // upward advection [MJ] int iter=0; double change=ALMOST_INF; double tolerance=1e-2; //[deg C] double T_guess,Ts_guess; double dE,dEs; - double E_guess =_aMres_last[p]; - double Es_guess =_aMsed_last[p]; + double E_guess =_aMres[p]; + double Es_guess =_aMsed[p]; while (change > tolerance) { - T_guess = ConvertVolumetricEnthalpyToTemperature(E_guess/V_e_new); + T_guess = ConvertVolumetricEnthalpyToTemperature(E_guess /V_e_new); Ts_guess= ConvertVolumetricEnthalpyToTemperature(Es_guess/V_h_new); A[0][0] = 1.0 + 0.5*tstep*(Q_new + Q_dn_new)/V_e_new; @@ -628,14 +555,102 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // E_guess +=dE; Es_guess +=dEs; iter++; - //if (iter>2){cout<<"iter: "<2){cout<<"iter: "<25){ + // Estimate from most recent known values and best guess for current temperatures + T_guess = ConvertVolumetricEnthalpyToTemperature(_aMres[p]/V_e_new); + Ts_guess= ConvertVolumetricEnthalpyToTemperature(_aMsed[p]/V_h_new); + dE = aMout_new[nSegments-1]*tstep; // inflow [MJ] + dE -= Q_new/V_e_new*_aMres[p]*tstep; // outflow [MJ] + dE -= Q_dn_new/V_e_new*_aMres[p]*tstep; // downward advection [MJ] + dE += Q_up_new/V_h_new*_aMsed[p]*tstep; // upward advection [MJ] + dE += _aMresRain[p]*tstep; // rainfall inputs [MJ] + LW = -STEFAN_BOLTZ*EMISS_WATER*pow(T_guess+ZERO_CELSIUS,4); + dE += A_new*(SW+LW_in+LW)*tstep; // net radiation [MJ] + dE -= A_new*(AET*DENSITY_WATER*LH_VAPOR)*tstep; // latent heat [MJ] + dE += A_new*hstar*(temp_air-T_guess)*tstep; // sensible heat exhange [MJ] + dE += kdiff*A_h_new*(Ts_guess-T_guess)*tstep; // diffusive heat exchange [MJ] + + dEs = Q_dn_new/V_e_new*_aMres[p]*tstep; // downward advection [MJ] + dEs+= Q_up_new/V_h_new*_aMsed[p]*tstep; // upward advection [MJ] + dEs+= kdiff*A_h_new*(T_guess-Ts_guess)*tstep; // diffusive heat exchange [MJ] + + E_guess = _aMres[p] + dE; + Es_guess = _aMsed[p] + dEs; + + break; + } } + string thisdate=tt.date_string; + string thishour=DecDaysToHours(tt.julian_day); + double stage = pRes->GetResStage (); + double T_new =ConvertVolumetricEnthalpyToTemperature(E_guess / V_e_new); + double Ts_new=ConvertVolumetricEnthalpyToTemperature(Es_guess / V_h_new ); + + /* + // Diagnostic data dump + cout << "Name ="<GetReservoirName()<SetOldDownwardFlowRate(Q_dn_new); + pRes->SetOldUpwardFlowRate(Q_up_new); delete [] B; delete [] R; @@ -1297,10 +1312,8 @@ void CEnthalpyModel::WriteOutputFileHeaders(const optStruct& Options) _LAKEOUT << name << " Eout [MJ/m2/d],"; _LAKEOUT << name << " Q_rain [MJ/m2/d],"; _LAKEOUT << name << " Q_sens [MJ/m2/d],"; - _LAKEOUT << name << " Q_cond [MJ/m2/d],"; + _LAKEOUT << name << " Q_conv [MJ/m2/d],"; _LAKEOUT << name << " Q_adv [MJ/m2/d],"; - _LAKEOUT << name << " Ri [0..1],"; - _LAKEOUT << name << " N2 [s-2],"; _LAKEOUT << name << " kdiff [MJ/m2/K/d],"; _LAKEOUT << name << " Q_lat [MJ/m2/d],"; _LAKEOUT << name << " Q_sw_in [MJ/m2/d],"; @@ -1334,10 +1347,10 @@ void CEnthalpyModel::WriteMinorOutput(const optStruct& Options,const time_struct // StreamReachEnergyBalances.csv //-------------------------------------------------------------------- - int p; + int p; double Q_sens,Q_lat,Q_GW,Q_sw_in,Q_lw_in,Q_lw_out,Q_fric,Q_rain,Q_adv,Tave; double Q_cond, Q_conv, Q_lateral; - double Ri, N2, kdiff; + double Ri, N2, kdiff; double Ein,Eout,mult; CSubBasin *pSB; @@ -1398,23 +1411,22 @@ void CEnthalpyModel::WriteMinorOutput(const optStruct& Options,const time_struct mult = 1.0 / HRUarea; - Ein = 0.5 * (_aMout_last[p] + _aMout[p][_pModel->GetSubBasin(p)->GetNumSegments() - 1]); - Eout = 0.5 * (_aMout_res [p] + _aMout_res_last[p]); + Ein = _aMout[p][_pModel->GetSubBasin(p)->GetNumSegments() - 1]; + Eout = _aMout_res [p]; double Qin=_pModel->GetSubBasin(p)->GetOutflowArray()[_pModel->GetSubBasin(p)->GetNumSegments()-1]*SEC_PER_DAY; //[m3/d] - GetEnergyLossesFromLake(p,Q_sens,Q_conv,Q_lat,Q_sw_in,Q_lw_in,Q_lw_out,Q_rain,Q_adv,Ri,N2,kdiff); + GetEnergyLossesFromLake(p,Q_sens,Q_conv,Q_lat,Q_sw_in,Q_lw_in,Q_lw_out,Q_rain,Q_adv,kdiff); double lakeTemp =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_e_new); double sedTemp =ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_h_new); double inflowTemp=ConvertVolumetricEnthalpyToTemperature( Ein / Qin ); double pctFroz =ConvertVolumetricEnthalpyToIceContent (_aMres[p] / V_e_new); -// double stage2_mixing = _stage - _mixing_depth; _LAKEOUT << mult * Ein << "," << mult * Eout << ","; _LAKEOUT << mult * Q_rain << "," << mult * Q_sens << ","; _LAKEOUT << mult * Q_conv << "," << mult * Q_adv << ","; - _LAKEOUT << Ri << "," << N2 << "," << kdiff << ","; + _LAKEOUT << kdiff << ","; _LAKEOUT << mult * Q_lat << ","; _LAKEOUT << mult * Q_sw_in << "," << mult * Q_lw_in << "," << mult * Q_lw_out << ","; _LAKEOUT << mult * _aMres[p] << "," << mult * _aMsed[p] << ","; diff --git a/src/EnergyTransport.h b/src/EnergyTransport.h index 411f02bb..b4565d9f 100644 --- a/src/EnergyTransport.h +++ b/src/EnergyTransport.h @@ -59,7 +59,7 @@ class CEnthalpyModel :public CConstituentModel double GetWaterTemperature (const double *state_vars, const int iWater) const; double GetEnergyLossesInTransit(const int p,double &Q_sens,double &Q_GW) const; - double GetEnergyLossesFromLake (const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out, double &Q_rain, double &Q_adv, double &Ri, double &N2, double &kdiff) const; + double GetEnergyLossesFromLake (const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out, double &Q_rain, double &Q_adv, double &kdiff) const; double GetEnergyLossesFromReach(const int p,double &Q_sens,double &Q_cond,double &Q_lat,double &Q_GW,double &Q_rad_in,double &Q_lw_in, double &Q_lw_out,double &Q_lateral, double &Q_fric, double &Tave) const; double GetOutflowIceFraction (const int p) const; diff --git a/src/Reservoir.cpp b/src/Reservoir.cpp index eb5443af..021e4338 100755 --- a/src/Reservoir.cpp +++ b/src/Reservoir.cpp @@ -24,6 +24,7 @@ void CReservoir::BaseConstructor(const string Name,const long SubID) _lakebed_thick=1.0; _lakebed_cond =0.0; _lake_convcoeff=2.0; + _lake_lyr_convcoeff=0.0; _stage =0.0; _stage_last=0.0; @@ -518,6 +519,11 @@ double CReservoir::GetLakebedConductivity() const { return _lakebed_cond; } // double CReservoir::GetLakeConvectionCoeff() const { return _lake_convcoeff; } +////////////////////////////////////////////////////////////////// +/// \returns lake between-layer thermal convection coefficient [MJ/m2/d/K] +// +double CReservoir::GetLakeLyrConvectionCoeff() const { return _lake_lyr_convcoeff; } + ////////////////////////////////////////////////////////////////// /// \returns current outflow rate [m3/s] // @@ -548,6 +554,17 @@ double CReservoir::GetIntegratedControlOutflow(const int i, const double& tstep { return 0.5*(_aQstruct[i]+_aQstruct_last[i])*(tstep*SEC_PER_DAY); //integrated } + +////////////////////////////////////////////////////////////////// +/// \returns previous downward flow rate between reservoir layers [m3/d] +// +double CReservoir::GetOldDownwardFlowRate() const { return _Q_dn_old; } + +////////////////////////////////////////////////////////////////// +/// \returns previous upward flow rate between reservoir layers [m3/d] +// +double CReservoir::GetOldUpwardFlowRate() const { return _Q_up_old; } + ////////////////////////////////////////////////////////////////// /// \returns previous storage [m3] // @@ -1031,6 +1048,12 @@ void CReservoir::SetLakeConvectionCoeff(const double& conv) { _lake_convcoeff=conv; } ////////////////////////////////////////////////////////////////// +/// \brief sets lake between-layer convection coeff +// +void CReservoir::SetLakeLyrConvectionCoeff(const double& conv) { + _lake_lyr_convcoeff=conv; +} +////////////////////////////////////////////////////////////////// /// \brief gets current constraint name /// \returns current constraint applied to estimate stage/flow // @@ -1511,6 +1534,24 @@ void CReservoir::SetMinStage(const double &min_z) _min_stage=min_z; } +////////////////////////////////////////////////////////////////// +/// \brief sets previous downward flow rate +/// \param Q_dn [in] downward flow rate between lake layers [m3/d] +// +void CReservoir::SetOldDownwardFlowRate(const double &Q_dn) +{ + _Q_dn_old=Q_dn; +} + +////////////////////////////////////////////////////////////////// +/// \brief sets previous upward flow rate +/// \param Q_up [in] upward flow rate between lake layers [m3/d] +// +void CReservoir::SetOldUpwardFlowRate(const double &Q_up) +{ + _Q_up_old=Q_up; +} + ////////////////////////////////////////////////////////////////// /// \brief Routes water through reservoir /// \param Qin_old [in] inflow at start of timestep diff --git a/src/Reservoir.h b/src/Reservoir.h index f4707c41..5c93606b 100755 --- a/src/Reservoir.h +++ b/src/Reservoir.h @@ -60,6 +60,7 @@ class CReservoir double _lakebed_thick; ///< lakebed thickness [m] double _lakebed_cond; ///< lakebed thermal conductivity [MJ/m/K/d] double _lake_convcoeff; ///< lake surface thermal convection coefficient [MJ/m2/d/K] + double _lake_lyr_convcoeff; ///< lake between-layer thermal convection coefficient [MJ/m2/d/K] const CHydroUnit *_pHRU; ///< (potentially zero-area) HRU used for Precip/ET calculation (or NULL for no ET) @@ -108,8 +109,8 @@ class CReservoir double _stage_last; ///< stage at beginning of current time step [m] double _Qout; ///< outflow corresponding to current stage [m3/s] double _Qout_last; ///< outflow at beginning of current time step [m3/s] - double _Q_dn_old; - double _Q_up_old; + double _Q_dn_old; ///< downward flow between layers at beginning of current time step [m3/d] + double _Q_up_old; ///< upward flow between layers at beginning of current time step [m3/d] double _MB_losses; ///< losses over current time step [m3] double _AET; ///< losses through AET only [m3] double _Precip; ///< gains through precipitation [m3] @@ -191,6 +192,8 @@ class CReservoir double GetOldOutflowRate () const; //[m3/s] double GetIntegratedOutflow (const double &tstep) const; //[m3] double GetControlOutflow (const int i) const; //[m3] + double GetOldUpwardFlowRate () const; //[m3/d] + double GetOldDownwardFlowRate () const; //[m3/d] double GetIntegratedControlOutflow(const int i, const double& tstep) const; //[m3] double GetReservoirLosses (const double &tstep) const; //[m3] double GetReservoirEvapLosses (const double &tstep) const; //[m3] @@ -205,6 +208,7 @@ class CReservoir double GetLakebedThickness () const; //[m] double GetLakebedConductivity () const; //[MJ/m/K/d] double GetLakeConvectionCoeff () const; //[MJ/m2/d/K] + double GetLakeLyrConvectionCoeff() const; //[MJ/m2/d/K] double GetDischargeFromStage (const double &stage, const int nn) const; //[m3/s] double GetStageDischargeDerivative(const double &stage,const int nn) const; //[m3/s/d] double GetMinStage (const int nn) const;//[m] @@ -250,6 +254,9 @@ class CReservoir void SetLakebedThickness (const double &thick); void SetLakebedConductivity (const double &cond); void SetLakeConvectionCoeff (const double &conv); + void SetLakeLyrConvectionCoeff(const double &conv); + void SetOldUpwardFlowRate (const double &Q_up); + void SetOldDownwardFlowRate (const double &Q_dn); void AddDemand (CDemand *pDemand); From c40065e9963ffd1251d5dc0c72c4c44cba35ba4e Mon Sep 17 00:00:00 2001 From: Markus Schnorbus Date: Tue, 11 Feb 2025 15:30:03 -0800 Subject: [PATCH 9/9] Added global parameter to set maximum value of convection/diffusion coefficient between lake layers, KDIFF_MAX. --- src/EnergyTransport.cpp | 77 +++++---------------------------------- src/GlobalParams.cpp | 15 +++++++- src/ParsePropertyFile.cpp | 3 ++ src/Properties.h | 1 + 4 files changed, 27 insertions(+), 69 deletions(-) diff --git a/src/EnergyTransport.cpp b/src/EnergyTransport.cpp index 9703b0b8..c14ce929 100755 --- a/src/EnergyTransport.cpp +++ b/src/EnergyTransport.cpp @@ -299,8 +299,8 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double Acorr=1.0; double SW(0), LW(0), LW_in(0), temp_air(0), AET(0); double hstar(0), ksed(0), Vsed=0.001; - double T_new =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_e_new); - double Ts_new=ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_h_new ); + double T_new =ConvertVolumetricEnthalpyToTemperature(_aMres[p] / V_e_new); + double Ts_new=ConvertVolumetricEnthalpyToTemperature(_aMsed[p] / V_h_new ); if(pHRU!=NULL) { //otherwise only simulate advective mixing+rain input Acorr =pHRU->GetArea()*M2_PER_KM2/A_new; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area @@ -369,6 +369,7 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // double A_h_old = pRes->GetOldMixingArea(); double A_h_avg = 0.5 * (A_h_new + A_h_old); double zm = pRes->GetMixingDepth(); + double kdiff_max = _pModel->GetGlobalParams()->GetParams()->lake_lyr_conv_coeff_max; double tstep=Options.timestep; @@ -470,8 +471,7 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // } //Final diffusion coefficient - kdiff = threshMin(md*(Ke+Ked+Km)*HCP_WATER, 5.0, 0.0); //[MJ/m2/K/d] - //kdiff = md*(Ke+Ked+Km)*HCP_WATER; //[MJ/m2/K/d] + kdiff = threshMin(md*(Ke+Ked+Km)*HCP_WATER, kdiff_max, 0.0); //[MJ/m2/K/d] //////////////////////////////// pRes->SetLakeLyrConvectionCoeff(kdiff); @@ -556,7 +556,7 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // Es_guess +=dEs; iter++; - if (iter>2){cout<<"iter: "<2){cout<<"iter: "<25){ @@ -585,67 +585,10 @@ void CEnthalpyModel::RouteMassInReservoir (const int p, // } } - string thisdate=tt.date_string; - string thishour=DecDaysToHours(tt.julian_day); - double stage = pRes->GetResStage (); - double T_new =ConvertVolumetricEnthalpyToTemperature(E_guess / V_e_new); - double Ts_new=ConvertVolumetricEnthalpyToTemperature(Es_guess / V_h_new ); - - /* - // Diagnostic data dump - cout << "Name ="<GetReservoirName()<25){cout<<"Name: "<GetReservoirName()<<" Date: "<