Skip to content

Commit

Permalink
updated bowen ratio estimation
Browse files Browse the repository at this point in the history
  • Loading branch information
Silwood authored and Silwood committed Oct 27, 2023
1 parent 42c56ee commit a9d1e2c
Show file tree
Hide file tree
Showing 10 changed files with 300 additions and 83 deletions.
6 changes: 6 additions & 0 deletions .vscode/settings.json
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
{
"files.associations": {
"cmath": "cpp"
},
"r.lsp.promptToInstall": false
}
4 changes: 2 additions & 2 deletions R/splash.point.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -366,9 +366,9 @@ soil_hydro<-function(sand, clay, OM, fgravel=0,bd=NA, ...) {
# 09. calc Air entry pressure [mm] from calibrated Saxton and Rawls (2006)
########################################################################################

moist_fvol33init<-0.278*sand+0.034*clay+0.022*OM-0.018*(sand*OM)-0.027*(clay*OM)-0.584*(sand*clay)+0.078
moist_fvol33init<-0.278*fsand+0.034*fclay+0.022*fOM-0.018*(fsand*fOM)-0.027*(fclay*fOM)-0.584*(fsand*fclay)+0.078
moist_fvol33<-moist_fvol33init+(0.636*moist_fvol33init-0.107)
bub_init<--21.6*sand-27.93*clay-81.97*moist_fvol33+71.12*(sand*moist_fvol33)+8.29*(clay*moist_fvol33)+14.05*(sand*clay)+27.16
bub_init<--21.6*fsand-27.93*fclay-81.97*moist_fvol33+71.12*(fsand*moist_fvol33)+8.29*(fclay*moist_fvol33)+14.05*(fsand*fclay)+27.16
bubbling_p<-bub_init+(0.02*bub_init^2-0.113*bub_init-0.7)
# 101.97162129779 converts from KPa to mmH2O
bubbling_p<-bubbling_p*-101.97162129779
Expand Down
5 changes: 5 additions & 0 deletions src/.vscode/settings.json
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
{
"files.associations": {
"cmath": "cpp"
}
}
81 changes: 48 additions & 33 deletions src/EVAP.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ void EVAP::calculate_daily_fluxes(double sw, int n, int y, double sw_in,
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 1. Calculate radiation values
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
solar.calculate_daily_fluxes(n, y, sw_in, tc, slop, asp,snow, nd);
solar.calculate_daily_fluxes(n, y, sw_in, tc, slop, asp,snow, nd, sw);
d_sr = solar.get_vals();
ru = d_sr.ru;
rv = d_sr.rv;
Expand Down Expand Up @@ -112,9 +112,11 @@ void EVAP::calculate_daily_fluxes(double sw, int n, int y, double sw_in,
lv = enthalpy_vap(tc);
pw = density_h2o(tc, patm);
g = psychro(tc, patm);
//econ = s/(lv*pw*(s + g));
// specific heat J/kg/K
//double Cp = specific_heat(tw);
econ = s/(lv*pw*(s + g));
// max evaporation from Yang & Roderick (2019) doi:10.1002/qj.3481
econ = s/(lv*pw*(s + 0.24*g));
//econ = s/(lv*pw*(s + 0.24*g));
visc = calc_viscosity_h2o(tw,patm);
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 3. Calculate daily condensation (wc), mm assume 10% of rnn_d (Jones, 2013)
Expand All @@ -124,7 +126,7 @@ void EVAP::calculate_daily_fluxes(double sw, int n, int y, double sw_in,
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 4. Estimate daily equilibrium evapotranspiration (eet_d), mm
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
eet_d = (1.0e3)*econ*rn_d;
eet_d = (1.0e3)*(s/(lv*pw*(s + 0.24*g)))*rn_d;

// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 5. Estimate daily potential evapotranspiration (pet_d), mm
Expand All @@ -143,15 +145,31 @@ void EVAP::calculate_daily_fluxes(double sw, int n, int y, double sw_in,
rx = (3.6e6)*econ;
////// calculate cw, assuming lim=1 at max demand (cos(h) =1)
// maximum instatntaneous demand (pet_max), mm/hr
double pet_max = rx*((rw*(ru+rv))- rnl) ;
pet_max = rx*((rw*(ru+rv))- rnl) ;
//assume cw = pet_max pet_max*0.55;
//sw *= pet_max*0.6;
sw *= pet_max*0.4;
//sw *= pet_max;
//sw *= 0.5;

// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 9. Estimate daily water supply, mm
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//// calc bowen ratio using eq. 17 from Gentine et al., (2011) doi:10.1175/2011JHM1261.1
double B_r = g/(sw*s);
//double B_r = (1/sw) + g/(sw*s) - 1.0;

double EF = 1/(B_r + 1.0);
///update sw
sw = pet_max*EF;

if (sw < 0.0 || isnan(sw)==1) {
sw = 0.0;
}
//else if (sw > pet_max){
// sw = pet_max;
//}
////////////////////////////////////////////////////////////////////////////////////////
//9.1. Soil matric potential to MPa
//double psi_m_mpa = psi_m * 0.00000980665;
//double psi_m_mpa = sw;
// //9.2. Transpiration potential to mol/m2
// double T_mol = pet_d / (18/pw);
// //9.3. Leaf water potential [atm]
Expand All @@ -169,21 +187,21 @@ void EVAP::calculate_daily_fluxes(double sw, int n, int y, double sw_in,
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/*
//9.1. Soil matric potential to MPa
double psi_m_mpa = psi_m * 0.00000980665;
double psi_m_mpa = sw;
// //9.2. Leaf water potential at critical leaf RWC [MPa]
// molar volume
double v_m = 18/pw;
double psi_l_c = ((0.082*(273.15+tw)*log(0.90))/(v_m))*0.101325;
// //9.4. Minimum Resistance soil to leaf asuming field capacity and Leaf water potential at critical leaf RWC [MPa]
double Rp = (-0.033-psi_l_c)/pet_d;
// //9.3. Leaf water potential assuming RWC at 98% [MPa]
double psi_l = ((0.082*(273.15+tw)*log(0.98))/(v_m))*0.101325;
//9.5. water supply mm/day
double sw = (psi_m_mpa-psi_l)/Rp;
double psi_l = ((0.082*(273.15+ts)*log(0.98))/(v_m))*0.101325;
//total resistance
// double Rp = (-1.0*psi_l)/pet_max;
// dimensioles Psi_s eq. 9.17 Campbell, 1998
double psi_ratio = psi_m_mpa/psi_l;
sw = (1.0 - (2.0 * psi_ratio/3.0))*pet_max;
sw *= 0.3;
*/

// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 7. Calculate the intersection hour angle (hi), degrees
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Expand All @@ -202,23 +220,16 @@ void EVAP::calculate_daily_fluxes(double sw, int n, int y, double sw_in,
// 8. Estimate daily snowmelt, mm and energy available for sublimation assume first snow melt then evaporate
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// estimate cold content [J/m^2] Hock, (2005) doi:10.1191/0309133305pp453ra [J/kg]
double Qsnow = 0.0;
if ((ts < 0.0)){
Qsnow = Global::Cpsnow * abs(ts) ;
}
// total energy required to raise the temperature and melt per kg [J/kg]
double Qmelt = Global::kfus + Qsnow;

//double Qsnow = 0.0;
//if ((ts < 0.0) && (elv > 1500)){
// Qsnow = Global::Cpsnow * pw * (snow/1000.0) * abs(ts) ;
//}
//snowmelt = min(snow,(rn_d/(pw*(Qmelt)))*1000.0);
// energy available for snowmelt
// if ((tc < 0.0) && (elv > 2000) && (snow > 0.0)){
// Qsnow = Global::Cpsnow * abs(tc) ;
// }
// available energy after heating the snow
//double AE = max(0.0,rn_d-Qsnow);
// calc snowmelt
//snowmelt = min(snow,(AE/(pw*Global::kfus))*1000.0);


// calc snowmelt

if (tc >= 3.0) {
snowmelt = min(snow,(rn_d/(pw*Global::kfus))*1000.0);
}else{
Expand All @@ -228,11 +239,13 @@ void EVAP::calculate_daily_fluxes(double sw, int n, int y, double sw_in,

// energy used in melting
double melt_enrg = (snowmelt/1000)*pw*Global::kfus;
// update available energy
//AE -= melt_enrg;
// energy available after snowmelt
double AE = rn_d - melt_enrg;
//calc evaporation from melted water (sublimation)
sublimation = min(snowmelt,(AE*econ)*1000.0);
// energy used in melting + evaporation
// energy used in heating + melting + evaporating snow
melt_enrg += ((sublimation/1000.0)/econ);

// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Expand All @@ -244,6 +257,7 @@ void EVAP::calculate_daily_fluxes(double sw, int n, int y, double sw_in,
aet_d *= (24.0/Global::PI);
//aet_d = sw*pet_d;
aet_d -= (melt_enrg*econ*1000.0);

if (aet_d < 0.0){
aet_d = 0.0;
}
Expand Down Expand Up @@ -517,6 +531,7 @@ etr EVAP::get_vals(){
d_etr.pw = pw;
d_etr.rn_d = rn_d;
d_etr.visc = visc;
d_etr.pet_max = pet_max;

return d_etr;
}
Expand Down
1 change: 1 addition & 0 deletions src/EVAP.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,7 @@ class EVAP {
double snowmelt; // daily snowmelt mm/d
double sublimation; // daily sublimation mm/d
double visc; // viscosity Pa s
double pet_max; // maximum (midday) instantaneous evapotranspiration mm/h
srad d_sr; // daily srad struct
etr d_etr; // daily etr struct

Expand Down
19 changes: 13 additions & 6 deletions src/SOLAR.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ SOLAR::SOLAR(double a, double b)
Class Function Definitions
/////////////////////////////////////////////////////////////////////
*/
void SOLAR::calculate_daily_fluxes(int n, int y, double sw_in, double tc, double slop, double asp,double snow, double nd){
void SOLAR::calculate_daily_fluxes(int n, int y, double sw_in, double tc, double slop, double asp,double snow, double nd, double sw){
/* ***********************************************************************
Name: SOLAR.calculate_daily_fluxes
Input: - int, day of year (n)
Expand All @@ -89,7 +89,9 @@ void SOLAR::calculate_daily_fluxes(int n, int y, double sw_in, double tc, double
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// @TODO: check for erroneous values
day = n;

if (sw > 1.0){
sw = 1.0;
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// 1. Calculate number of days in year (kN), days
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Expand Down Expand Up @@ -209,9 +211,12 @@ void SOLAR::calculate_daily_fluxes(int n, int y, double sw_in, double tc, double
// lm with modis and snotel data, max snow albedo as 0.78
//double sfc = 0.5935119/(1.0 + exp((-243.1452366 - snow)/173.6612577));
double sfc =snow/(140.0+snow);
//calc albedo as function of soil moisture,
double alb_v = Global::alb_sw-0.17*sw;
//albedo, assuming max snow albedo as 0.85

double alb = Global::alb_sw *(1.0-sfc)+sfc*max_alb_snw;
//double alb = Global::alb_sw *(1.0-sfc)+sfc*max_alb_snw;
double alb = alb_v *(1.0-sfc)+sfc*max_alb_snw;
if ((sw_in == 0.0) || (hs == 0.0)){
rw = (1.0 - alb)*tau*dr*Global::Gsc;
} else {
Expand Down Expand Up @@ -252,9 +257,11 @@ void SOLAR::calculate_daily_fluxes(int n, int y, double sw_in, double tc, double
// 14. Calculate mean surface temperature (C) Yang & Roderick (2019) doi:10.1002/qj.3481
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//double LW_in =
ts = tc + 2.517*exp(2.38*tau) + 0.03466 *abs(lat) ;
//ts = tc + 2.517*exp(2.38*tau) + 0.03466 *abs(lat) ;
// empirical fitting Fluxnet
//ts = 1.0691*tc + 2.7302*tau - 0.035*abs(lat) + 2.334;
ts = 1.0691*tc + 2.7302*tau - 0.035*abs(lat) + 2.334;
//max air temeperature
//ts = tc/dcos(hs/2);


}
Expand Down Expand Up @@ -380,7 +387,7 @@ srad SOLAR::get_vals(){
dsr.hn = hn; // cross-over hour angle, degrees
dsr.rn_d = rn_d; // daytime net radiation, J/m^2
dsr.rnn_d = rnn_d; // nighttime net radiation, J/m^2
dsr.ts = ts; // surface temperature, J/m^2
dsr.ts = ts; // surface temperature, C

return dsr;
}
Expand Down
2 changes: 1 addition & 1 deletion src/SOLAR.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ class SOLAR {
SOLAR(double a, double b);

// General purpose functions:
void calculate_daily_fluxes(int n, int y, double sw_in, double tc, double slop, double asp, double snow, double nd);
void calculate_daily_fluxes(int n, int y, double sw_in, double tc, double slop, double asp, double snow, double nd, double sw);
srad get_vals();
void display();
};
Expand Down
Loading

0 comments on commit a9d1e2c

Please sign in to comment.