From a9d1e2c747048006c50b32d20a99505977eeb0bb Mon Sep 17 00:00:00 2001 From: Silwood Date: Fri, 27 Oct 2023 14:07:02 +0100 Subject: [PATCH] updated bowen ratio estimation --- .vscode/settings.json | 6 + R/splash.point.R | 4 +- src/.vscode/settings.json | 5 + src/EVAP.cpp | 81 +++++++----- src/EVAP.h | 1 + src/SOLAR.cpp | 19 ++- src/SOLAR.h | 2 +- src/SPLASH.cpp | 257 ++++++++++++++++++++++++++++++++------ src/etr.h | 1 + src/global.cpp | 7 +- 10 files changed, 300 insertions(+), 83 deletions(-) create mode 100644 .vscode/settings.json create mode 100644 src/.vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..a5db189 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,6 @@ +{ + "files.associations": { + "cmath": "cpp" + }, + "r.lsp.promptToInstall": false +} \ No newline at end of file diff --git a/R/splash.point.R b/R/splash.point.R index cc989cc..e9aad00 100644 --- a/R/splash.point.R +++ b/R/splash.point.R @@ -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 diff --git a/src/.vscode/settings.json b/src/.vscode/settings.json new file mode 100644 index 0000000..1be854f --- /dev/null +++ b/src/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "files.associations": { + "cmath": "cpp" + } +} \ No newline at end of file diff --git a/src/EVAP.cpp b/src/EVAP.cpp index 62cb2b7..d6de0c2 100644 --- a/src/EVAP.cpp +++ b/src/EVAP.cpp @@ -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; @@ -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) @@ -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 @@ -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] @@ -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 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -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{ @@ -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); // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -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; } @@ -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; } diff --git a/src/EVAP.h b/src/EVAP.h index 36108bb..92ca8ae 100644 --- a/src/EVAP.h +++ b/src/EVAP.h @@ -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 diff --git a/src/SOLAR.cpp b/src/SOLAR.cpp index e1f6a91..3f85097 100644 --- a/src/SOLAR.cpp +++ b/src/SOLAR.cpp @@ -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) @@ -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 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -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 { @@ -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); } @@ -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; } diff --git a/src/SOLAR.h b/src/SOLAR.h index 38ed7c7..71448fa 100644 --- a/src/SOLAR.h +++ b/src/SOLAR.h @@ -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(); }; diff --git a/src/SPLASH.cpp b/src/SPLASH.cpp index 8923985..6c99c33 100644 --- a/src/SPLASH.cpp +++ b/src/SPLASH.cpp @@ -184,7 +184,7 @@ void SPLASH::quick_run(int n, int y, double wn, double sw_in, double tc, //bubbling pressure/capillarity fringe (mm) double bub_press = soil_info[6]; //residual water content, test as WP? - //double RES = WP/2; + //double RES = WP; double RES = soil_info[7]; //upslope area double Au = soil_info[8]; @@ -260,6 +260,21 @@ void SPLASH::quick_run(int n, int y, double wn, double sw_in, double tc, // #################################################################################################################### // 03. calculate supply rate (sw) // #################################################################################################################### + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // 1. Calculate evaporative supply rate (sw), mm/h Original + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + double max_sw = 1.0/ pow((1 + pow(AI,6.8681)),0.07956432); + // // org splash + double sw = ((wn-RES)/(Wmax-RES)); + //double sw = ((wn)/(Wmax)); + if (sw < 0.0 || isnan(sw)==1) { + sw = 0.0; + } else if (sw > 1.0){ + sw = 1.0; + } + /* double sw = 0.0; @@ -397,15 +412,24 @@ void SPLASH::quick_run(int n, int y, double wn, double sw_in, double tc, //if (sw < 0.0 || isnan(sw)==1) { // sw = 0.0; //} + */ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // Water supply calc as difference of soil and water potentials mm // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + /* - - evap.calculate_daily_fluxes(0.0, n, y, sw_in, tc, slop, asp,snow,0.0); + evap.calculate_daily_fluxes(1.0, n, y, sw_in, tc, slop, asp,snow,nd); etr dn_dumm = evap.get_vals(); - double pet_d = dn_dumm.pet; + //double pet_d = dn_dumm.pet; double pw = dn_dumm.pw; + visc = dn_dumm.visc; + // 5.1.1. correct Ksat for fuid viscosity + double int_perm = Ksat/Global::fluidity; + double Ksat_visc = int_perm*((pw*Global::G)/visc)*3.6; + // soil conductivity unsaturated [mm/h] + double Kunsat = Ksat_visc * pow((theta_i/theta_s),(3.0+(2.0/lambda))); + // maximum (midday) potential evapotranspiration mm/h + double pet_max = dn_dumm.pet_max; //9.1. Soil matric potential to MPa double psi_m_mpa = psi_m * 0.00000980665; @@ -417,26 +441,56 @@ void SPLASH::quick_run(int n, int y, double wn, double sw_in, double tc, //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+tc)*log(0.98))/(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.02-psi_l)/pet_d; + double psi_l = ((0.082*(278.15+tc)*log(0.95))/(v_m))*0.101325; + //double psi_l = -1*pow(-4*psi_m_mpa,0.5); + //9.4. Minimum Resistance soil to leaf asuming SATURATION and Leaf water potential at critical leaf RWC [MPa] + // assuming path water 2 m, then water head is [MPa] + //double h_path = pw*Global::G*2.0*1E-6; + //////////////////////////////////////////////////////// + ///calc potential conductivity + // total max conductivity [mm/h]= root_leaf K + soilK + //double K_T = pet_max * h_path/(-1.0*psi_l); + //vertical hidraulic gradient at saturation + //double up_hg = (0.033/h_path) - 1.0; + //double K_T = pet_max/up_hg ; + // maximum root to leaf conductivity [mm/h] + //double K_RL = (K_T*Ksat_visc)/(Ksat_visc-K_T); + // correct K_RL for viscosity + //double K_RL_cte = K_RL * visc; + //////////////////////////////////////////////////////// + ///calc real conductivity k_plant=0.5^(((-6:0)/-4)^3) + //double KRL_h = K_RL* pow(0.5,pow(psi_l/-4,3)); + //double KT_real = (K_RL*Kunsat)/(K_RL+Kunsat); + //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; - //9.5. water supply mm/day - sw = (psi_m_mpa-psi_l)/Rp; - //sw = (1.0 - (2.0 * psi_ratio/3.0))* pet_d; + //9.5. water supply mm/h + //double sw = (psi_m_mpa-psi_l)/Rp; + //9.5. water supply mm/h + //double sw = KT_real * (psi_m_mpa-psi_l)/h_path; + //double sw = KT_real * ( ((psi_m_mpa-psi_l)/h_path)- 1.0); + ///use sw to input in et module + //double sw = psi_m_mpa; + + double sw = (1.0 - (1/3)*(psi_ratio))*pet_max; + if (sw < 0.0 || isnan(sw)==1) { sw = 0.0; } + //double max_sw = 0.7/ pow((1 + pow(AI,6.8681)),0.07956432); + //sw *= max_sw; */ - - + //////////////////// package + + /* double max_sw = 1.0/ pow((1 + pow(AI,6.8681)),0.07956432); double sw = 0.0; - + double upt = 0.0 ; // when the soil depth exeeds 2m: if (depth>=6.0){ @@ -445,18 +499,53 @@ void SPLASH::quick_run(int n, int y, double wn, double sw_in, double tc, } else{ // bedrock < 2 m - sw = ((wn-RES)/(Wmax-RES)); + //upt = ((wn-RES)/(Wmax-RES)); + upt = ((wn-WP)/(SAT-WP)); + sw = 1.0 - pow(1.0 + 0.5*upt,-1.0/lambda); } - + //sw *= max_sw; + sw *= Global::Cw; if (sw < 0.0 || isnan(sw)==1) { sw = 0.0; - } else if (sw > 1.0){ - sw = 1.0; - } + } + + //else if (sw > max_sw){ + // sw = max_sw; + //} - sw *= max_sw; - + */ + + + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // new theory combined + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + /* + //parameters aridity ceiling function + double m_AI = 6.8681; + double n_AI = 0.07956432; + // double max_sw = 1.0/ pow((1 + pow(AI,m_AI)),n_AI); + ///shape limiting function + double g_AI = pow( 1 - pow(1 + pow(AI,m_AI),-1*n_AI), -1.0*lambda)- 1.0; + // available water content in fraction + double awc = (wn)/(Wmax); + /// fix bounds + if (awc < 0.0 || isnan(awc)==1) { + awc = 0.0; + } + //else if (awc > 1.0){ + // awc = 1.0; + //} + // calc fraction of the water uptake + double sw = 1.0 - pow(1 + g_AI*awc, -1.0/lambda); + + if (sw < 0.0 || isnan(sw)==1) { + sw = 0.0; + } + //else if (sw > max_sw){ + // sw = max_sw; + //} + */ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -868,7 +957,7 @@ void SPLASH::run_one_day(int n, int y, double wn, double sw_in, double tc, //bubbling pressure/capillarity fringe (mm) double bub_press = soil_info[6]; //residual water content, test as WP? - //double RES = WP/2; + //double RES = WP; double RES = soil_info[7]; //upslope area double Au = soil_info[8]; @@ -942,6 +1031,22 @@ void SPLASH::run_one_day(int n, int y, double wn, double sw_in, double tc, // #################################################################################################################### // 03. calculate supply rate (sw) // #################################################################################################################### + + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // 1. Calculate evaporative supply rate (sw), mm/h Original + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + double max_sw = 1.0/ pow((1 + pow(AI,6.8681)),0.07956432); + // // org splash + double sw = ((wn-RES)/(Wmax-RES)); + //double sw = ((wn)/(Wmax)); + if (sw < 0.0 || isnan(sw)==1) { + sw = 0.0; + } else if (sw > 1.0){ + sw = 1.0; + } + /* double sw = 0.0; // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -978,15 +1083,23 @@ void SPLASH::run_one_day(int n, int y, double wn, double sw_in, double tc, } else if (sw > Global::Cw){ sw = Global::Cw; } - + */ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // Water supply calc as difference of soil and water potentials mm // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - evap.calculate_daily_fluxes(0.0, n, y, sw_in, tc, slop, asp,snow,0.0); + /* + evap.calculate_daily_fluxes(1.0, n, y, sw_in, tc, slop, asp,snow,0.0); etr dn_dumm = evap.get_vals(); double pet_d = dn_dumm.pet; double pw = dn_dumm.pw; + visc = dn_dumm.visc; + // 5.1.1. correct Ksat for fuid viscosity + double int_perm = Ksat/Global::fluidity; + double Ksat_visc = int_perm*((pw*Global::G)/visc)*3.6; + // soil conductivity unsaturated [mm/h] + double Kunsat = Ksat_visc * pow((theta_i/theta_s),(3.0+(2.0/lambda))); + // maximum (midday) potential evapotranspiration mm/h + double pet_max = dn_dumm.pet_max; //9.1. Soil matric potential to MPa double psi_m_mpa = psi_m * 0.00000980665; @@ -998,25 +1111,57 @@ void SPLASH::run_one_day(int n, int y, double wn, double sw_in, double tc, //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+tc)*log(0.98))/(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.02-psi_l)/pet_d; + double psi_l = ((0.082*(278.15+tc)*log(0.95))/(v_m))*0.101325; + //double psi_l = -1*pow(-4*psi_m_mpa,0.5); + //9.4. Minimum Resistance soil to leaf asuming SATURATION and Leaf water potential at critical leaf RWC [MPa] + // assuming path water 2 m, then water head is [MPa] + //double h_path = pw*Global::G*2.0*1E-6; + //////////////////////////////////////////////////////// + ///calc potential conductivity + // total max conductivity [mm/h]= root_leaf K + soilK + //double K_T = pet_max * h_path/(-1.0*psi_l); + //vertical hidraulic gradient at saturation + //double up_hg = (0.033/h_path) - 1.0; + //double K_T = pet_max/up_hg ; + // maximum root to leaf conductivity [mm/h] + //double K_RL = (K_T*Ksat_visc)/(Ksat_visc-K_T); + // correct K_RL for viscosity + //double K_RL_cte = K_RL * visc; + //////////////////////////////////////////////////////// + ///calc real conductivity k_plant=0.5^(((-6:0)/-4)^3) + //double KRL_h = K_RL* pow(0.5,pow(psi_l/-4,3)); + //double KT_real = (K_RL*Kunsat)/(K_RL+Kunsat); + //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; - //9.5. water supply mm/day - sw = (psi_m_mpa-psi_l)/Rp; - //sw = (1.0 - (2.0 * psi_ratio/3.0))*pet_d; + //9.5. water supply mm/h + //double sw = (psi_m_mpa-psi_l)/Rp; + //9.5. water supply mm/h + //double sw = KT_real * (psi_m_mpa-psi_l)/h_path; + //double sw = KT_real * ( ((psi_m_mpa-psi_l)/h_path)- 1.0); + ///use sw to input in et module + //double sw = psi_m_mpa; + + double sw = (1.0 - (1/3)*(psi_ratio))*pet_max; + if (sw < 0.0 || isnan(sw)==1) { sw = 0.0; } + //double max_sw = 0.7/ pow((1 + pow(AI,6.8681)),0.07956432); + //sw *= max_sw; */ + //////////////////// package + /* + double max_sw = 1.0/ pow((1 + pow(AI,6.8681)),0.07956432); double sw = 0.0; - + double upt = 0.0 ; // when the soil depth exeeds 2m: if (depth>=6.0){ @@ -1024,18 +1169,54 @@ void SPLASH::run_one_day(int n, int y, double wn, double sw_in, double tc, sw = ((w_z-RES_z)/(Wmax-RES_z)); } else{ // bedrock < 2 m - - sw = ((wn-RES)/(Wmax-RES)); + //upt = ((wn-RES)/(Wmax-RES)); + upt = ((wn-WP)/(SAT-WP)); + sw = 1.0 - pow(1.0 + 0.5*upt,-1.0/lambda); } + //sw *= max_sw; + sw *= Global::Cw; if (sw < 0.0 || isnan(sw)==1) { sw = 0.0; - } else if (sw > 1.0){ - sw = 1.0; - } + } + //else if (sw > max_sw){ + // sw = max_sw; + // } + */ + - sw *= max_sw; + + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // new theory combined + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + /* + //parameters aridity ceiling function + double m_AI = 6.8681; + double n_AI = 0.07956432; + //double max_sw = 1.0/ pow((1 + pow(AI,m_AI)),n_AI); + ///shape limiting function + double g_AI = pow( 1 - pow(1 + pow(AI,m_AI),-1*n_AI), -1.0*lambda)- 1.0; + // available water content in fraction + double awc = (wn)/(Wmax); + /// fix bounds + if (awc < 0.0 || isnan(awc)==1) { + awc = 0.0; + } + //else if (awc > 1.0){ + // awc = 1.0; + // } + // calc fraction of the water uptake + double sw = 1.0 - pow(1 + g_AI*awc, -1.0/lambda); + + if (sw < 0.0 || isnan(sw)==1) { + sw = 0.0; + } + + //else if (sw > max_sw){ + // sw = max_sw; + //} + */ // #################################################################################################################### // 04. Snowpack and energy Balances // #################################################################################################################### @@ -1438,7 +1619,7 @@ List SPLASH::spin_up(int n, int y, vector &sw_in, vector &tair, //bubbling pressure/capillarity fringe (mm) double bub_press = soil_info[6]; //residual water content, test as WP? - //double RES = WP/2; + //double RES = WP; double RES = soil_info[7]; double theta_s = SAT/(depth *1000.0); diff --git a/src/etr.h b/src/etr.h index 4e31ab2..8aae90c 100644 --- a/src/etr.h +++ b/src/etr.h @@ -52,5 +52,6 @@ struct etr { double pw; // density of water, kg/m^3 double rn_d; // daytime net radioation J/m^2 double visc; // Viscosity Pa s + double pet_max; // maximum (midday) instantaneous evapotranspiration mm/h }; #endif diff --git a/src/global.cpp b/src/global.cpp index 1478231..caa9820 100644 --- a/src/global.cpp +++ b/src/global.cpp @@ -49,7 +49,8 @@ namespace Global { //extern const double A(170); // (Monteith & Unsworth, 1990) extern const double A(91.86328); // (Sandoval et al., 2021) - extern const double alb_sw(0.17); // (Federer, 1968) + extern const double alb_sw(0.30); // (Stephens, 2015) + //extern const double alb_sw(0.17); // (Federer, 1968) extern const double alb_vis(0.03); // (Sellers, 1985) //extern const double b(0.20); // (Linacre, 1968) extern const double b(0.2012435); // (Sandoval et al., 2021) @@ -67,8 +68,8 @@ namespace Global { extern const double To(288.15); // (Berberan-Santos et al., 1997) extern const double w(0.26); // (Priestley & Taylor, 1972) //extern const double w(-0.3); // (Priestley & Taylor, 1972) - extern const double Cw(1.0); // (Federer, 1982) - //extern const double Cw(0.3); // (Federer, 1982) + //extern const double Cw(1.05); // (Federer, 1982) + extern const double Cw(0.45); // (Federer, 1982) extern const double Wm(180.0); // (Cramer & Prentice, 1988) extern const double PI(3.141592653589793); // pi extern const double pir = (PI/180.0); // degrees to radians