From 416502c127dd538ef32765f2abb8b0f27f855beb Mon Sep 17 00:00:00 2001 From: ZSong Date: Mon, 5 Aug 2024 15:03:22 +0200 Subject: [PATCH 01/17] Change R_depth, soil depth, soil layers thickness for CH-HTC site --- src/+io/getModelSettings.m | 6 +-- src/Dtrmn_Z.m | 108 ++++++++++++++++++++++--------------- 2 files changed, 69 insertions(+), 45 deletions(-) diff --git a/src/+io/getModelSettings.m b/src/+io/getModelSettings.m index b39cc6a6..c36f6c66 100644 --- a/src/+io/getModelSettings.m +++ b/src/+io/getModelSettings.m @@ -2,7 +2,7 @@ %{ %} - ModelSettings.R_depth = 350; + ModelSettings.R_depth = 30; % Indicator denotes the index of soil type for choosing soil physical parameters ModelSettings.J = 1; @@ -51,7 +51,7 @@ ModelSettings.rroot = 1.5 * 1e-3; ModelSettings.SFCC = 1; - ModelSettings.Tot_Depth = 500; % Unit is cm. it should be usually bigger than 0.5m. Otherwise, + ModelSettings.Tot_Depth = 100; % Unit is cm. it should be usually bigger than 0.5m. Otherwise, ModelSettings.Eqlspace = 0; % Indicator for deciding is the space step equal or not; % the DeltZ would be reset in 50cm by hand; ModelSettings.NS = 1; % Number of soil types; @@ -61,7 +61,7 @@ ModelSettings.KT = 0; % Number of time steps; % Determination of NL, the number of elments - ModelSettings.NL = 100; + ModelSettings.NL = 29; if ~ModelSettings.Eqlspace [DeltZ, DeltZ_R, NL, ML] = Dtrmn_Z(ModelSettings.NL, ModelSettings.Tot_Depth); else diff --git a/src/Dtrmn_Z.m b/src/Dtrmn_Z.m index cb0a3145..76ace63d 100644 --- a/src/Dtrmn_Z.m +++ b/src/Dtrmn_Z.m @@ -3,51 +3,75 @@ The determination of the element length %} Elmn_Lnth = 0; + sitename = 'CH-HTC'; + if strcmp(sitename,'CH-HTC') + for ML=1:6 + DeltZ_R(ML)=1;%4 + end + % DeltZ_R(6)=1;%4 + for ML=7:13 + DeltZ_R(ML)=2;%5 + end + for ML=14:21 + DeltZ_R(ML)=2.5;%5 + end + for ML=22:25 + DeltZ_R(ML)=5; + end + for ML=26:29 + DeltZ_R(ML)=10; + end + NL= ML; + for ML=1:NL + MML=NL-ML+1; + DeltZ(ML)=DeltZ_R(MML); + end + else + for ML = 1:3 + DeltZ_R(ML) = 1; % 4 + end + DeltZ_R(4) = 2; % 4 + for ML = 5:14 + DeltZ_R(ML) = 2; % 5 + end + for ML = 15:18 + DeltZ_R(ML) = 2.5; % 5 + end + for ML = 19:23 + DeltZ_R(ML) = 5; + end + for ML = 24:31 + DeltZ_R(ML) = 10; + end + for ML = 32:40 + DeltZ_R(ML) = 10; + end + for ML = 41:42 + DeltZ_R(ML) = 15; % 5 + end + % Sum of element lengths and compared to the total lenght, so that judge + % can be made to determine the length of rest elements. - for ML = 1:3 - DeltZ_R(ML) = 1; % 4 - end - DeltZ_R(4) = 2; % 4 - for ML = 5:14 - DeltZ_R(ML) = 2; % 5 - end - for ML = 15:18 - DeltZ_R(ML) = 2.5; % 5 - end - for ML = 19:23 - DeltZ_R(ML) = 5; - end - for ML = 24:31 - DeltZ_R(ML) = 10; - end - for ML = 32:40 - DeltZ_R(ML) = 10; - end - for ML = 41:42 - DeltZ_R(ML) = 15; % 5 - end - % Sum of element lengths and compared to the total lenght, so that judge - % can be made to determine the length of rest elements. - - for ML = 1:42 - Elmn_Lnth = Elmn_Lnth + DeltZ_R(ML); - end + for ML = 1:42 + Elmn_Lnth = Elmn_Lnth + DeltZ_R(ML); + end - % If the total sum of element lenth is over the predefined depth, stop the - % for loop, make the ML, at which the element lenth sumtion is over defined - % depth, to be new NL. - DeltZ = []; - for ML = 43:NL - DeltZ_R(ML) = 20; - Elmn_Lnth = Elmn_Lnth + DeltZ_R(ML); - if Elmn_Lnth >= Tot_Depth - DeltZ_R(ML) = Tot_Depth - Elmn_Lnth + DeltZ_R(ML); - NL = ML; - for ML = 1:NL - MML = NL - ML + 1; - DeltZ(ML) = DeltZ_R(MML); + % If the total sum of element lenth is over the predefined depth, stop the + % for loop, make the ML, at which the element lenth sumtion is over defined + % depth, to be new NL. + DeltZ = []; + for ML = 43:NL + DeltZ_R(ML) = 20; + Elmn_Lnth = Elmn_Lnth + DeltZ_R(ML); + if Elmn_Lnth >= Tot_Depth + DeltZ_R(ML) = Tot_Depth - Elmn_Lnth + DeltZ_R(ML); + NL = ML; + for ML = 1:NL + MML = NL - ML + 1; + DeltZ(ML) = DeltZ_R(MML); + end + return end - return end end end From 68adc700d04ca48cd5db25b078be184ca84c705f Mon Sep 17 00:00:00 2001 From: ZSong Date: Sun, 11 Aug 2024 22:42:57 +0200 Subject: [PATCH 02/17] Medlyn stomatal conductance --- src/+io/loadTimeSeries.m | 6 +- src/+io/select_input.m | 4 + src/+parameters/loadParameters.m | 2 +- src/+parameters/setTempParameters.m | 30 +++++- src/STEMMUS_SCOPE.m | 5 +- src/biochemical.m | 147 ++++++++++++++++++++++------ src/ebal.m | 10 ++ 7 files changed, 168 insertions(+), 36 deletions(-) diff --git a/src/+io/loadTimeSeries.m b/src/+io/loadTimeSeries.m index 96329f87..0155f438 100644 --- a/src/+io/loadTimeSeries.m +++ b/src/+io/loadTimeSeries.m @@ -122,7 +122,7 @@ ScopeParameters.SMC = load(fullfile(path_input, SMC_file)); end - %% 7. Set Leaf Vcmo, Tparam, m, Type, Rdparam, and leafwidth as time-dependent + %% 7. Set Leaf Vcmo, Tparam, m, Type, Rdparam,leafwidth, g1Med and g0Med as time-dependent % parameters. The timeseries of landcover, along with the lookup table values % (lc...) are used to generate the timeseries. landcovers = unique(landcoverClass, 'stable'); @@ -134,9 +134,11 @@ ScopeParameters.Type(timeIndex, 1) = ScopeParameters.lcType(landcoverIndex, 1); ScopeParameters.Rdparam(timeIndex, 1) = ScopeParameters.lcRdparam(landcoverIndex, 1); ScopeParameters.leafwidth(timeIndex, 1) = ScopeParameters.lcleafwidth(landcoverIndex, 1); + % ScopeParameters.g1Med(timeIndex, 1) = ScopeParameters.lcg1Med(landcoverIndex, 1); + % ScopeParameters.g0Med(timeIndex, 1) = ScopeParameters.lcg0Med(landcoverIndex, 1); end % Remove the intermediate variables that were required to generate the timeseries: - ScopeParameters = rmfield(ScopeParameters, {'lcVcmo', 'lcTparam', 'lcm', 'lcType', 'lcRdparam', 'lcleafwidth'}); + ScopeParameters = rmfield(ScopeParameters, {'lcVcmo', 'lcTparam', 'lcm', 'lcType', 'lcRdparam', 'lcleafwidth', 'lcg1Med', 'lcg0Med'}); if ~isempty(Cab_file) Cabtable = load(fullfile(path_input, Cab_file)); diff --git a/src/+io/select_input.m b/src/+io/select_input.m index b91810e4..bc8e5adc 100644 --- a/src/+io/select_input.m +++ b/src/+io/select_input.m @@ -30,6 +30,8 @@ leafbio.Tparam = ScopeParameters.Tparam(digitsVector(14), :); % this is correct (: instead of 14) fqe = ScopeParameters.fqe(digitsVector(15)); leafbio.Rdparam = ScopeParameters.Rdparam(digitsVector(13)); + leafbio.g1Med = ScopeParameters.g1Med(digitsVector(65)); + leafbio.g0Med = ScopeParameters.g0Med(digitsVector(66)); leafbio.rho_thermal = ScopeParameters.rho_thermal(digitsVector(7)); leafbio.tau_thermal = ScopeParameters.tau_thermal(digitsVector(8)); @@ -39,6 +41,8 @@ leafbio.kNPQs = ScopeParameters.kNPQs(digitsVector(57)); leafbio.qLs = ScopeParameters.qLs(digitsVector(58)); leafbio.stressfactor = ScopeParameters.stressfactor(digitsVector(59)); + leafbio.g1Med = ScopeParameters.g1Med; + leafbio.g0Med = ScopeParameters.g0Med; canopy.LAI = ScopeParameters.LAI(digitsVector(22)); canopy.hc = ScopeParameters.hc(digitsVector(23)); diff --git a/src/+parameters/loadParameters.m b/src/+parameters/loadParameters.m index 499d74ab..d7e3efdf 100644 --- a/src/+parameters/loadParameters.m +++ b/src/+parameters/loadParameters.m @@ -12,7 +12,7 @@ 'Ca', 'Oa', 'rb', 'Cd', 'CR', 'CD1', 'Psicor', 'CSSOIL', 'rbs', 'rwc', ... 'startDOY', 'endDOY', 'LAT', 'LON', 'timezn', 'tts', 'tto', 'psi', 'SMC', ... 'Tyear', 'beta', 'kNPQs', 'qLs', 'stressfactor', 'Cant', 'BSMBrightness', ... - 'BSMlat', 'BSMlon', 'BallBerry0' + 'BSMlat', 'BSMlon', 'BallBerry0', 'g1Med', 'g0Med' }; % create an empty structure with field names of ScopeParametersNames diff --git a/src/+parameters/setTempParameters.m b/src/+parameters/setTempParameters.m index a72cc22d..f4dad17e 100644 --- a/src/+parameters/setTempParameters.m +++ b/src/+parameters/setTempParameters.m @@ -6,7 +6,7 @@ % where landcoverClass is an array like {"forest", "forest", "forest", "shrubland", ...} landcovers = unique(landcoverClass, 'stable'); for landcoverIndex = 1:length(landcovers) - [Vcmo, Tparam, m, Type, Rdparam, leafwidth] = landcover_variables(landcovers(landcoverIndex), siteName); + [Vcmo, Tparam, m, g1Med, g0Med, Type, Rdparam, leafwidth] = landcover_variables(landcovers(landcoverIndex), siteName); % The lc... parameters are intermediate values. The timeseries are generated in % io.loadTimeSeries.m. ScopeParameters.lcVcmo(landcoverIndex, 1) = Vcmo; @@ -15,41 +15,53 @@ ScopeParameters.lcType(landcoverIndex, 1) = Type; ScopeParameters.lcRdparam(landcoverIndex, 1) = Rdparam; ScopeParameters.lcleafwidth(landcoverIndex, 1) = leafwidth; + ScopeParameters.lcg1Med(landcoverIndex, 1) = g1Med; + ScopeParameters.lcg0Med(landcoverIndex, 1) = g0Med; end end %% lookup table for Vcmo, Tparam, m, Type, Rdparam, and leafwidth -function [Vcmo, Tparam, m, Type, Rdparam, leafwidth] = landcover_variables(landCoverType, siteName) +function [Vcmo, Tparam, m, g1Med, g0Med, Type, Rdparam, leafwidth] = landcover_variables(landCoverType, siteName) Rdparam = [0.015]; % NOTE: Rdparam IS MISSING FOR MOST DEFINTIONS: if startsWith(landCoverType, 'Permanent Wetlands') Tparam = [0.2 0.3 288 313 328]; % These are five parameters specifying the temperature response. Vcmo = [120]; % Vcmax, maximum carboxylation capacity (at optimum temperature) m = [9]; % Ball-Berry stomatal conductance parameter + g1Med = [14.43]; % Medlyn stomatal conductance parameter: g1 + g0Med = [0.001]; % Medlyn stomatal conductance parameter: g0 Type = [0]; % Photochemical pathway: 0=C3, 1=C4 leafwidth = [0.05]; % leaf width elseif startsWith(landCoverType, 'Evergreen Broadleaf') Tparam = [0.2 0.3 283 311 328]; Vcmo = [80]; m = [9]; + g1Med = [9]; + g0Med = [0.001]; Type = [0]; leafwidth = [0.05]; elseif startsWith(landCoverType, 'Deciduous Broadleaf') Tparam = [0.2 0.3 283 311 328]; Vcmo = [80]; m = [9]; + g1Med = [6.99]; + g0Med = [-0.044]; Type = [0]; leafwidth = [0.05]; elseif startsWith(landCoverType, 'Mixed Forests') Tparam = [0.2 0.3 281 307 328]; Vcmo = [80]; m = [9]; + g1Med = [5.36]; + g0Med = [0.024]; Type = [0]; leafwidth = [0.04]; elseif startsWith(landCoverType, 'Evergreen Needleleaf') Tparam = [0.2 0.3 278 303 328]; Vcmo = [80]; m = [9]; + g1Med = [1.66]; + g0Med = [0.033]; Type = [0]; leafwidth = [0.01]; elseif startsWith(landCoverType, 'Croplands') @@ -67,28 +79,38 @@ Type = [0]; leafwidth = [0.03]; end + g1Med = [9]; + g0Med = [0.001]; elseif startsWith(landCoverType, 'Open Shrublands') Tparam = [0.2 0.3 288 313 328]; Vcmo = [120]; m = [9]; + g1Med = [5.95]; + g0Med = [0.028]; Type = [0]; leafwidth = [0.05]; elseif startsWith(landCoverType, 'Closed Shrublands') Tparam = [0.2 0.3 288 313 328]; Vcmo = [80]; m = [9]; + g1Med = [5.95]; + g0Med = [0.028]; Type = [0]; leafwidth = [0.05]; elseif startsWith(landCoverType, 'Savannas') Tparam = [0.2 0.3 278 313 328]; Vcmo = [120]; m = [9]; + g1Med = [11.23]; + g0Med = [-0.004]; Type = [0]; leafwidth = [0.05]; elseif startsWith(landCoverType, 'Woody Savannas') Tparam = [0.2 0.3 278 313 328]; Vcmo = [120]; m = [9]; + g1Med = [11.23]; + g0Med = [-0.004]; Type = [0]; leafwidth = [0.03]; elseif startsWith(landCoverType, 'Grassland') @@ -104,12 +126,16 @@ Type = [0]; end leafwidth = [0.02]; + g1Med = [11.23]; + g0Med = [-0.004]; else Tparam = [0.2 0.3 288 313 328]; Vcmo = [80]; m = [9]; Type = [0]; leafwidth = [0.05]; + g1Med = [9]; + g0Med = [0.001]; warning('landCover name unknown, "%s" is not recognized. ', landCoverType); end end diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 69838e47..4cdd3d60 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -1,4 +1,4 @@ -%% STEMMUS-SCOPE.m (script) + %% STEMMUS-SCOPE.m (script) % STEMMUS-SCOPE is a model for Integrated modeling of canopy photosynthesis, fluorescence, % and the transfer of energy, mass, and momentum in the soil-plant-atmosphere continuum % Copyright (C) 2021 Yunfei Wang, Lianyu Yu, Yijian Zeng, Christiaan Van der Tol, Bob Su @@ -110,6 +110,7 @@ % Create a structure holding Scope parameters useXLSX = 1; % set it to 1 or 0, the current stemmus-scope does not support 0 [ScopeParameters, options] = parameters.loadParameters(options, useXLSX, X, F, N); + options.gsMethod = 2; % 1 for BallBerry's method; 2 for Medlyn's method % Define the location information ScopeParameters.LAT = SiteProperties.latitude; % latitude @@ -125,6 +126,8 @@ % Input T parameters for different vegetation type [ScopeParameters] = parameters.setTempParameters(ScopeParameters, SiteProperties.sitename, SiteProperties.landcoverClass); + ScopeParameters.g1Med = ScopeParameters.lcg1Med(1,1); + ScopeParameters.g0Med = ScopeParameters.lcg0Med(1,1); %% 5. Declare paths path_input = InputPath; % path of all inputs diff --git a/src/biochemical.m b/src/biochemical.m index 33e3aaa5..a5f79463 100644 --- a/src/biochemical.m +++ b/src/biochemical.m @@ -136,6 +136,11 @@ eb = biochem_in.eb; O = biochem_in.O; p = biochem_in.p; + ei = biochem_in.ei; + + gsMethod = biochem_in.gsMethod; + g1Med = biochem_in.g1Med; + g0Med = biochem_in.g0Med; % physiological Type = biochem_in.Type; @@ -318,47 +323,85 @@ %% calculation of Ci (internal CO2 concentration) RH = min(1, eb ./ equations.satvap(T - 273.15)); % jak: don't allow "supersaturated" air! (esp. on T curves) + VPD_l2b = max(0.0001, ei-eb); warnings = []; fcount = 0; % the number of times we called computeA() - if ~isempty(Ci_input) - Ci = Ci_input; % in units of bar. - if any(Ci_input > 1) - % assume Ci_input is in units of ppm. Convert to bar - Ci = Ci_input .* ppm2bar; - end - A = computeA(Ci); - elseif all(BallBerry0 == 0) - % b = 0: no need to iterate: - Ci = BallBerry(Cs, RH, [], BallBerrySlope, BallBerry0, minCi); - A = computeA(Ci); + switch gsMethod + case 1 % BallBerry's stomatal conductance + if ~isempty(Ci_input) + Ci = Ci_input; % in units of bar. + if any(Ci_input > 1) + % assume Ci_input is in units of ppm. Convert to bar + Ci = Ci_input .* ppm2bar; + end + A = computeA(Ci); + + elseif all(BallBerry0 == 0) + % b = 0: no need to iterate: + Ci = BallBerry(Cs, RH, [], BallBerrySlope, BallBerry0, minCi); + A = computeA(Ci); + + else + % compute Ci using iteration (JAK) + % it would be nice to use a built-in root-seeking function but fzero requires scalar inputs and outputs, + % Here I use a fully vectorized method based on Brent's method (like fzero) with some optimizations. + tol = 1e-6; % 0.1 ppm more-or-less + % Setting the "corner" argument to Gamma may be useful for low Ci cases, but not very useful for atmospheric CO2, so it's ignored. + % (fn, x0, corner, tolerance) + Ci = equations.fixedp_brent_ari(@(x) Ci_next(x, Cs, RH, minCi), Cs, [], tol); % [] in place of Gamma: it didn't make much difference + % NOTE: A is computed in Ci_next on the final returned Ci. fixedp_brent_ari() guarantees that it was done on the returned values. + % A = computeA(Ci); + end - else - % compute Ci using iteration (JAK) - % it would be nice to use a built-in root-seeking function but fzero requires scalar inputs and outputs, - % Here I use a fully vectorized method based on Brent's method (like fzero) with some optimizations. - tol = 1e-6; % 0.1 ppm more-or-less - % Setting the "corner" argument to Gamma may be useful for low Ci cases, but not very useful for atmospheric CO2, so it's ignored. - % (fn, x0, corner, tolerance) - Ci = equations.fixedp_brent_ari(@(x) Ci_next(x, Cs, RH, minCi), Cs, [], tol); % [] in place of Gamma: it didn't make much difference - % NOTE: A is computed in Ci_next on the final returned Ci. fixedp_brent_ari() guarantees that it was done on the returned values. - % A = computeA(Ci); + case 2 % Medlyn's stomatal conductance + if ~isempty(Ci_input) + Ci = Ci_input; % in units of bar. + if any(Ci_input > 1) + % assume Ci_input is in units of ppm. Convert to bar + Ci = Ci_input .* ppm2bar; + end + A = computeA(Ci); + + elseif all(g0Med == 0) + % b = 0: no need to iterate: + Ci = Medlyn(Cs, A, VPD_l2b, g1Med, g0Med, minCi); + A = computeA(Ci); + + else + % compute Ci using iteration (JAK) + % it would be nice to use a built-in root-seeking function but fzero requires scalar inputs and outputs, + % Here I use a fully vectorized method based on Brent's method (like fzero) with some optimizations. + tol = 1e-7; % 0.1 ppm more-or-less + % Setting the "corner" argument to Gamma may be useful for low Ci cases, but not very useful for atmospheric CO2, so it's ignored. + % (fn, x0, corner, tolerance) + Ci = equations.fixedp_brent_ari(@(x) Ci_next_ME(x, Cs, VPD_l2b, minCi), Cs, [], tol); % [] in place of Gamma: it didn't make much difference + %NOTE: A is computed in Ci_next on the final returned Ci. fixedp_brent_ari() guarantees that it was done on the returned values. + %A = computeA(Ci); + end end %% Test-function for iteration % (note that it assigns A in the function's context.) % As with the next section, this code can be read as if the function body executed at this point. % (if iteration was used). In other words, A is assigned at this point in the file (when iterating). - function [err, Ci_out] = Ci_next(Ci_in, Cs, RH, minCi) + function [err, Ci_out] = Ci_next_BB(Ci_in, Cs, RH, minCi) % compute the difference between "guessed" Ci (Ci_in) and Ci computed using BB after computing A - A = computeA(Ci_in); - A_bar = A .* ppm2bar; - Ci_out = BallBerry(Cs, RH, A_bar, BallBerrySlope, BallBerry0, minCi); % [Ci_out, gs] - + A = computeA(Ci_in); % A: ppm + A_bar = A .* ppm2bar; % A: bar + Ci_out = BallBerry(Cs, RH, A_bar, BallBerrySlope, BallBerry0, minCi); %[Ci_out, gs] Ci: bar + + err = Ci_out - Ci_in; % f(x) - x + end + function [err, Ci_out] = Ci_next_ME(Ci_in, Cs, VPD_l2b, minCi) + % compute the difference between "guessed" Ci (Ci_in) and Ci computed using BB after computing A + A = computeA(Ci_in); % A: ppm + A_bar = A .* ppm2bar; % A: bar + Ci_out = Medlyn(Cs, A_bar, VPD_l2b,g1Med, g0Med, minCi); %[Ci_out, gs] Ci: bar + err = Ci_out - Ci_in; % f(x) - x end - %% Compute Assimilation. % Note: even though computeA() is written as a separate function, % the code is, in fact, executed exactly this point in the file (i.e. between the previous if clause and the next section @@ -547,7 +590,7 @@ Ci = Ci_input; gs = []; if ~isempty(A) && nargout > 1 - gs = gsFun(Cs, RH, A, BallBerrySlope, BallBerry0); + gs = gsFun_BB(Cs, RH, A, BallBerrySlope, BallBerry0); end elseif all(BallBerry0 == 0) || isempty(A) % EXPLANATION: *at equilibrium* CO2_in = CO2_out => A = gs(Cs - Ci) [1] @@ -564,13 +607,13 @@ % note: the original B-B units are A: umol/m2/s, ci ppm (umol/mol), RH (unitless) % Cs input was ppm but was multiplied by ppm2bar above, so multiply A by ppm2bar to put them on the same scale. % don't let gs go below its minimum value (i.e. when A goes negative) - gs = gsFun(Cs, RH, A, BallBerrySlope, BallBerry0); + gs = gsFun_BB(Cs, RH, A, BallBerrySlope, BallBerry0); Ci = max(minCi .* Cs, Cs - 1.6 * A ./ gs); end end % function -function gs = gsFun(Cs, RH, A, BallBerrySlope, BallBerry0) +function gs = gsFun_BB(Cs, RH, A, BallBerrySlope, BallBerry0) % add in a bit just to avoid div zero. 1 ppm = 1e-6 (note since A < 0 if Cs ==0, it gives a small gs rather than maximal gs gs = max(BallBerry0, BallBerrySlope .* A .* RH ./ (Cs + 1e-9) + BallBerry0); % clean it up: @@ -578,6 +621,50 @@ gs(isnan(Cs)) = NaN; % max(NaN, X) = X (MATLAB 2013b) so fix it here end +%% Medlyn gs +function [Ci, gs] = Medlyn(Cs, A, VPD_l2b, g1Med, g0Med, minCi, Ci_input) + % Cs : CO2 at leaf surface + % A : Net assimilation in 'same units of CO2 as Cs'/m2/s + % VPD_l2b: vapour pressure deficit from leaf to leaf boundary [mbar or hPa] + % BallBerrySlope, BallBerry0: Medlynslope, Medlyn0 + % minCi : minimum Ci as a fraction of Cs (in case RH is very low?) + % Ci_input : will calculate gs if A is specified. + if nargin > 6 && ~isempty(Ci_input) + % Ci is given: try and compute gs + Ci = Ci_input; + gs = []; + if ~isempty(A) && nargout > 1 + gs = gsFun_ME(Cs, A, VPD_l2b, g1Med, g0Med); + end + elseif all(g0Med == 0) || isempty(A) + % EXPLANATION: *at equilibrium* CO2_in = CO2_out => A = gs(Cs - Ci) [1] + % so Ci = Cs - A/gs (at equilibrium) [2] + % Medlyn suggest: gs = g0+(1+g1/VPD^0.5)* A/Cs ( see. Medlyn 2010 GCB) + % % if g0 = 0 we can rearrange B-B for the second term in [2]: + + % % Ci = Cs - Cs/(VPD^0.5/(m+VPD^0.5)) = Cs ( 1- 1/(VPD^0.5/(m+VPD^0.5) [ the 1.6 converts from CO2- to H2O-diffusion ] + Ci = max(minCi .* Cs, Cs.*(1-(1.6.*VPD_l2b.^0.5./(g1Med + VPD_l2b.^0.5)))); + gs = []; + else + % % if b > 0 Ci = Cs( 1 - 1/(m RH + b Cs/A) ) + % % if we use Leuning 1990, Ci = Cs - (Cs - Gamma)/(m RH + b(Cs - Gamma)/A) [see def of Gamma, above] + % % note: the original B-B units are A: umol/m2/s, ci ppm (umol/mol), RH (unitless) + % % Cs input was ppm but was multiplied by ppm2bar above, so multiply A by ppm2bar to put them on the same scale. + % % don't let gs go below its minimum value (i.e. when A goes negative) + gs = gsFun_ME(Cs, A, VPD_l2b, g1Med, g0Med); + Ci = max(minCi .* Cs, Cs - 1.6 * A./gs) ; + end + +end % function + +function gs = gsFun_ME(Cs, A, VPD_l2b, g1Med, g0Med) + % add in a bit just to avoid div zero. 1 ppm = 1e-6 (note since A < 0 if Cs ==0, it gives a small gs rather than maximal gs + gs = max(g0Med, 1.6.*(1+(g1Med./ VPD_l2b .^ 0.5)) .* A ./ (Cs+1e-9) + g0Med ); + + %gs( Cs == 0 ) = would need to be max gs here; % eliminate infinities + gs( isnan(Cs) ) = NaN; % max(NaN, X) = X (MATLAB 2013b) so fix it here +end + %% Fluorescence model function [eta, qE, qQ, fs, fo, fm, fo0, fm0, Kn] = Fluorescencemodel(ps, x, Kp, Kf, Kd, Knparams) % note: x isn't strictly needed as an input parameter but it avoids code-duplication (of po0) and it's inherent risks. diff --git a/src/ebal.m b/src/ebal.m index d7af30c8..e29eff79 100644 --- a/src/ebal.m +++ b/src/ebal.m @@ -171,6 +171,8 @@ resistance of leaves (or biochemical_MD12: alternative) end LAI = canopy.LAI; + eih = equations.satvap(Tch); + eiu = equations.satvap(Tcu); PSI = 0; [bbx] = Max_Rootdepth(InitialValues.bbx); @@ -235,8 +237,14 @@ resistance of leaves (or biochemical_MD12: alternative) biochem_in.p = p; biochem_in.m = leafbio.m; biochem_in.BallBerry0 = leafbio.BallBerry0; + biochem_in.g1Med = leafbio.g1Med; + biochem_in.g0Med = leafbio.g0Med; + biochem_in.gsMethod = options.gsMethod; biochem_in.O = meteo.Oa; biochem_in.Rdparam = leafbio.Rdparam; + biochem_in.g1Med = leafbio.g1Med; + biochem_in.g0Med = leafbio.g0Med; + if options.Fluorescence_model == 2 % specific for the v.Caemmerer-Magnani model biochem_in.Tyear = leafbio.Tyear; @@ -257,6 +265,7 @@ resistance of leaves (or biochemical_MD12: alternative) biochem_in.Vcmo = fVh .* leafbio.Vcmo; biochem_in.Cs = Cch; biochem_in.Q = rad.Pnh_Cab * 1E6; + biochem_in.ei = eih; if options.Fluorescence_model == 2 biochem_out = biochemical_MD12(biochem_in); @@ -279,6 +288,7 @@ resistance of leaves (or biochemical_MD12: alternative) biochem_in.Vcmo = fVu .* leafbio.Vcmo; biochem_in.Cs = Ccu; biochem_in.Q = rad.Pnu_Cab * 1E6; + biochem_in.ei = eiu; if options.Fluorescence_model == 2 biochem_out = biochemical_MD12(biochem_in); From 74c664746e7deb0d89cc102a1d372f743af9915d Mon Sep 17 00:00:00 2001 From: ZSong Date: Tue, 13 Aug 2024 15:50:21 +0200 Subject: [PATCH 03/17] add multiple gs and PHS options --- src/STEMMUS_SCOPE.m | 10 ++++++++-- src/setScenario.m | 46 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 2 deletions(-) create mode 100644 src/setScenario.m diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 4cdd3d60..d9738c0e 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -40,6 +40,9 @@ % Read the configPath file. Due to using MATLAB compiler, we cannot use run(CFG) disp (['Reading config from ', CFG]); [InputPath, OutputPath, InitialConditionPath] = io.read_config(CFG); + gsOption = 'Medlyn'; % 1 for BallBerry's method; 2 for Medlyn's method + phsOption = 2; % 1 for CLM5, 2 for ED2, 3 for PHS + % Prepare forcing and soil data [SiteProperties, SoilProperties, TimeProperties] = io.prepareInputData(InputPath); @@ -68,6 +71,9 @@ RWUs = []; RWUg = []; + %% Set senarios + [biochemical, gsMethod, phwsfMethod] = setScenario(gsOption, phsOption); + %% 1. define Constants Constants = io.define_constants(); @@ -110,7 +116,7 @@ % Create a structure holding Scope parameters useXLSX = 1; % set it to 1 or 0, the current stemmus-scope does not support 0 [ScopeParameters, options] = parameters.loadParameters(options, useXLSX, X, F, N); - options.gsMethod = 2; % 1 for BallBerry's method; 2 for Medlyn's method + options.gsMethod = gsMethod; % 1 for BallBerry's method; 2 for Medlyn's method % Define the location information ScopeParameters.LAT = SiteProperties.latitude; % latitude @@ -352,7 +358,7 @@ % Will do one timestep in "update mode", and run until the end if in "full run" mode. while KT < endTime - KT = KT + 1; % Counting Number of timesteps + KT = KT + 1 % Counting Number of timesteps if KT > 1 && Delt_t > (TEND - TIME) Delt_t = TEND - TIME; % If Delt_t is changed due to excessive change of state variables, the judgement of the last time step is excuted. end diff --git a/src/setScenario.m b/src/setScenario.m new file mode 100644 index 00000000..25531713 --- /dev/null +++ b/src/setScenario.m @@ -0,0 +1,46 @@ +function [biochemical, gsMethod, phwsfMethod] = setScenario(gsOption, phsOption) +%{ + This function is used to set the stomatal conductance scheme and plant + hydraulic pathway option. + + input: + gsOption: a string indicate which stomatal conductance will be used + phsOption: a number indicate whether use plant hydraulic pathway + + output: + biochemical: a string, indicate which function will be used + gsMethod: a number indicate the stomatal conductance scheme: + 1 for BallBerry + 2 for Medlyn + phsOption: a number indicate if open the plant hydraulic pathway. + 1 for CLM5, + 2 for ED2, + 3 for PHS + +%} + + %% Set scenario + biochemical = @biochemical; + if strcmp(gsOption, 'BallBerry') + gsMethod = 1; + elseif strcmp(gsOption, 'Medlyn') + gsMethod = 2; + end + %% PHS setting + if phsOption % if phsOption is true, PHS open. Set plant water stress method. + switch phsOption + case 1 + phwsfMethod = 'CLM5'; + case 2 + phwsfMethod = 'ED2'; + case 3 + phwsfMethod = 'PHS'; + otherwise + fprintf('Please set phsOption\n'); + end + else % if phsOption is false, PHS close. + phwsfMethod = 'WSF'; + fprintf('PHS close, use soil water stress factor\n') + end + +end \ No newline at end of file From 09c949bc4a333aeaf7b170bb0c7d71abcf97e7f2 Mon Sep 17 00:00:00 2001 From: ZSong Date: Tue, 13 Aug 2024 17:13:48 +0200 Subject: [PATCH 04/17] define plant hydraulic constants --- src/+io/define_plant_constants.m | 255 +++++++++++++++++++++++++++++++ src/STEMMUS_SCOPE.m | 14 +- src/calRootProperties.m | 81 ++++++++++ 3 files changed, 348 insertions(+), 2 deletions(-) create mode 100644 src/+io/define_plant_constants.m create mode 100644 src/calRootProperties.m diff --git a/src/+io/define_plant_constants.m b/src/+io/define_plant_constants.m new file mode 100644 index 00000000..1172fcc4 --- /dev/null +++ b/src/+io/define_plant_constants.m @@ -0,0 +1,255 @@ +function ParaPlant = define_plant_constants(SiteProperties, phwsfMethod) +% define parameters used in plant hydraulic pathway +% Input: +% SiteProperties +% phwsfMethod: define phwsf method + +% Output: +% ParaPlant: A structure contains all of parameters used in plant hydraulic pathway. + +% Reference: +% Parameters cite from Li,Hongmei_forests_2021 +% Functions cite from Kennedy_JAMES_2019 + + + %% ----------------------- root growth ---------------------------- + ParaPlant.C2B = 2; %(g biomass / gC) + ParaPlant.B2C = 0.488; % ratio of biomass to carbon [gC / g biomass] + ParaPlant.rootRadius = 1.5*1e-3; % root radius [m] (2.9e-4 in CLM5), (0.5-6e-3 m)in STEMMUS_SCOPE_GMD + ParaPlant.rootDensity = 250*1000; % Root density Jackson et al., 1997 [gDM / m^3] + pa2m = 1/9810; % Pressure = rho g h = 1000kg m-3 * 9.81 m s-2 *1 m = 9810 Pa. 1m = 9810Pa + ParaPlant.s2l = 0.1; % The ratio of stem area index to leaf area index: SAI = s2l * LAI; + + %% ----------------------- define IGBP vegetation type ---------------------- + igbpVegLong = SiteProperties.landcoverClass; + ParaPlant.phwsfMethod = phwsfMethod; + switch phwsfMethod + case 'CLM5' + %% ----------------------- Plant hydraulics ----------------------- + if strcmp(igbpVegLong(1:18)', 'Permanent Wetlands') + ParaPlant.p50Leaf = -260; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -260; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -260; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 3.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule in CLM + ParaPlant.shapeFactorStem = 3.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 3.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =2e-8;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =2e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =2e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + + elseif strcmp(igbpVegLong(1:19)', 'Evergreen Broadleaf') + ParaPlant.p50Leaf = -260; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -175; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -175; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 1.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 1.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 1.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =6e-9;%1.28e-7;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =4e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =4e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + + elseif strcmp(igbpVegLong(1:19)', 'Deciduous Broadleaf') + ParaPlant.p50Leaf = -270; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -270; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -270; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 3.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 3.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 3.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =2e-8;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =2e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =2e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + + elseif strcmp(igbpVegLong(1:13)', 'Mixed Forests') + ParaPlant.p50Leaf = -260; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -260; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -260; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 3.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 3.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 3.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =2e-8;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =2e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =2e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + + elseif strcmp(igbpVegLong(1:20)', 'Evergreen Needleleaf') + ParaPlant.p50Leaf = -465; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -465; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -465; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 3.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 3.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 3.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =2e-8;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =2e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =2e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + + elseif strcmp(igbpVegLong(1:9)', 'Croplands') + ParaPlant.p50Leaf = -340; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -340; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -340; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 3.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 3.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 3.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =2e-8;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =2e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =2e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + + elseif strcmp(igbpVegLong(1:15)', 'Open Shrublands') + ParaPlant.p50Leaf = -260; % parameter of plant hydraulic pathway [m]. + ParaPlant.p50Stem = -260; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -260; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 3.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 3.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 3.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =2e-8;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =2e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =2e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + + elseif strcmp(igbpVegLong(1:17)', 'Closed Shrublands') + ParaPlant.p50Leaf = -260; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -260; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -260; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 3.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 3.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 3.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =2e-8;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =2e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =2e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + + elseif strcmp(igbpVegLong(1:8)', 'Savannas') + ParaPlant.p50Leaf = -340; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -340; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -340; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 3.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 3.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 3.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =2e-8;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =2e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =2e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + + elseif strcmp(igbpVegLong(1:14)', 'Woody Savannas') + ParaPlant.p50Leaf = -340; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -340; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -340; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 3.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 3.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 3.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =2e-8;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =2e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =2e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + + elseif strcmp(igbpVegLong(1:9)', 'Grassland') + ParaPlant.p50Leaf = -340; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -340; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -340; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 3.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 3.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 3.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =2e-8;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =2e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =2e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + else + warning("Don't find IGBP vegetation type! \n") + end + case 'ED2' + if strcmp(igbpVegLong(1:19)', 'Evergreen Broadleaf') + ParaPlant.p50Leaf = -67.67; %-175; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -342; %-60; %-175; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -30; %-175; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 1.98; %2.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 1.98; %0.8; %2.95; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 0.8; %2.95; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =2e-9;%1.28e-7;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =4e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =4e-8;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + else + + ParaPlant.p50Leaf = -67.67; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -342.51; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -300; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 1.98; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 1.98; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 1.98; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =1.28e-7;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =3.88e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =1.77e-10;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + end + case 'PHS' + if strcmp(igbpVegLong(1:19)', 'Evergreen Broadleaf') + ParaPlant.p50Leaf = -67.67; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Stem = -342.51; % parameter of plant hydraulic pathway [m] + ParaPlant.p50Root = -300; % parameter of plant hydraulic pathway [m] + + ParaPlant.shapeFactorLeaf = 1.98; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM + ParaPlant.shapeFactorStem = 1.98; % parameter of plant hydraulic pathway [unitless] + ParaPlant.shapeFactorRoot = 1.98; % parameter of plant hydraulic pathway [unitless] + + ParaPlant.Krootmax =1.28e-7;%2e-9; % root conductivity [m s-1] + ParaPlant.Kstemmax =3.88e-8;%2e-8; % stem conductivity [m s-1] + ParaPlant.Kleafmax =1.77e-10;%2e-7; % maximum leaf conductance [s-1] + ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + end + otherwise + fprintf('phwsf method need to be defined.') + end + + + +% % ParaPlant.psi50_sunleaf = ; % parameter of plant hydraulic pathway [MPa] +% % ParaPlant.psi50_shdleaf = ; % parameter of plant hydraulic pathway [MPa] +% % ParaPlant.ck_sunleaf = ; % parameter of plant hydraulic pathway [MPa] +% % ParaPlant.ck_shdleaf = ; % parameter of plant hydraulic pathway [MPa] +% +% ParaPlant.p50Leaf = -260;%50.97;%-1.75; %0.5; % unit MPa 5.0968e-11; % 0.5./1e6*pa2m; % parameter of plant hydraulic pathway [m] 0.5MPa refer to Kennedy 2019. +% ParaPlant.p50Stem = -260;%-1.75; % unit MPa -1.7839e-10; % -1.75./1e6*pa2m; % parameter of plant hydraulic pathway [m] -1.75MPa refer to Kennedy 2019. +% ParaPlant.p50Root = -260;%-1.75; % unit MPa -1.7839e-10; % -1.75./1e6*pa2m; % parameter of plant hydraulic pathway [m] -1.75MPa refer to Kennedy 2019. +% +% ParaPlant.ckLeaf = 3.95; % parameter of plant hydraulic pathway, 2.95 is the default vaule inCLM +% ParaPlant.ckStem = 3.95; % parameter of plant hydraulic pathway [unitless] +% ParaPlant.ckRoot = 3.95; % parameter of plant hydraulic pathway [unitless] +% +% ParaPlant.Krootmax =2e-8;%2e-9; % root conductivity [m s-1] +% ParaPlant.Kstemmax =2e-8;%2e-8; % stem conductivity [m s-1] +% ParaPlant.Kleafmax =2e-8;%2e-7; % maximum leaf conductance [s-1] +% ParaPlant.rootLateralLength = 0.25; % average coarse root length [m] + + +end + diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index d9738c0e..07af247b 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -16,7 +16,7 @@ % % You should have received a copy of the GNU General Public License % along with this program. If not, see . - +clear all; clc; % Load in required Octave packages if STEMMUS-SCOPE is being run in Octave: if exist('OCTAVE_VERSION', 'builtin') ~= 0 disp('Loading Octave packages...'); @@ -76,10 +76,20 @@ %% 1. define Constants Constants = io.define_constants(); + ParaPlant = io.define_plant_constants(SiteProperties, phwsfMethod); RTB = 1000; % initial root total biomass (g m-2) + if strcmp(SiteProperties.sitename, 'CH-HTC') + RTB = 300; + end % Rl used in ebal - [Rl, Ztot] = Initial_root_biomass(RTB, ModelSettings.DeltZ_R, ModelSettings.rroot, ModelSettings.ML, SiteProperties.landcoverClass(1)); + % [Rl, Ztot] = Initial_root_biomass(RTB, ModelSettings.DeltZ_R, ModelSettings.rroot, ModelSettings.ML, SiteProperties.landcoverClass(1)); + numSoilLayer = ModelSettings.ML; + soilThickness = ModelSettings.DeltZ_R'; + [RootProperties,soilDepth] = calRootProperties(SiteProperties, ParaPlant, numSoilLayer, soilThickness, RTB); + Rl = RootProperties.lengthDensity; + Ztot = soilDepth; + %% 2. simulation options path_input = InputPath; % path of all inputs diff --git a/src/calRootProperties.m b/src/calRootProperties.m new file mode 100644 index 00000000..0d274da3 --- /dev/null +++ b/src/calRootProperties.m @@ -0,0 +1,81 @@ +function [RootProperties, soilDepth] = calRootProperties(SiteProperties, ParaPlant, numSoilLayer, soilThickness, RTB) +% This function is used to calculate root properties. +% Input: +% SiteProperties: A structure contains site properties. In this function, vegetation type (igbpVegLong) is needed. +% ParaPlant : A structure contains plant parameters +% numSoilLayer : The number of soil layers +% soilTickness : soil thickness [cm] +% RTB : root to biomass [g m-2]. This parameter is defined in +% Constants.m with the value of 1000. +% +% Output: +% RootProperties: A structure contains root properties. + + + %% ============= define plant parameters ============== + B2C = ParaPlant.B2C; %[g C / g biomass] + rootRadius = ParaPlant.rootRadius; %[m] + rootDensity = ParaPlant.rootDensity; %[gDM / m^3] + igbpVegLong = SiteProperties.landcoverClass; + + %% Jackson et al. (1996), it is refered to CLM5.0 document. + if strcmp(igbpVegLong(1:18), 'Permanent Wetlands') + beta = 0.993; + elseif strcmp(igbpVegLong(1:19)', 'Evergreen Broadleaf') + beta = 0.993; + elseif strcmp(igbpVegLong(1:19)', 'Deciduous Broadleaf') + beta = 0.993; + elseif strcmp(igbpVegLong(1:13)', 'Mixed Forests') + beta = 0.993; + elseif strcmp(igbpVegLong(1:20)', 'Evergreen Needleleaf') + beta = 0.993; + elseif strcmp(igbpVegLong(1:9)', 'Croplands') + beta = 0.943; + elseif strcmp(igbpVegLong(1:15)', 'Open Shrublands') + beta = 0.966; + elseif strcmp(igbpVegLong(1:17)', 'Closed Shrublands') + beta = 0.966; + elseif strcmp(igbpVegLong(1:8)', 'Savannas') + beta = 0.943; + elseif strcmp(igbpVegLong(1:14)', 'Woody Savannas') + beta = 0.943; + elseif strcmp(igbpVegLong(1:9)', 'Grassland') + beta = 0.943; + else + beta = 0.943; + warning('IGBP vegetation name unknown, "%s" is not recognized. Falling back to default value for beta', IGBP_veg_long) + end + + %% ======= calculate the depth of each soil layer to land surface ===== + % change the unit from cm to m + soilDepth=cumsum(soilThickness,1); % [unit: cm] + + %% ============= calculate root distribution ================== + for i=1:numSoilLayer + if i==1 + rootFrac(i)=(1-beta.^(soilDepth(i)/2)); % root fraction in each soil layer + else + rootFrac(i)=beta.^(soilDepth(i-1)-(soilThickness(i-1)/2))-beta.^(soilDepth(i-1)+(soilThickness(i)/2)); + end + end + rootFrac = rootFrac'; + + %% ================ root cross area =========================== + rootCrossArea = pi .* rootRadius^2; % [m^2] + + %% =============== root length density ====================== + % RTB = 1000, gC/m2 or g biomass/ m2? I think the unit is gC/m2. + Rltot = RTB/B2C/rootDensity/rootCrossArea; % root length index [m root / m^2 PFT] + rootLengthDensity = (Rltot .* rootFrac) ./ (soilThickness./100); % [m/m3]; 100 is scale factor from cm to m. + + %% ==================== root spacing ======================= + % root spacing is same with CLM5 based on eq.11.12 in CLM5 tech notes. + rootSpac = sqrt(1./(pi.* rootLengthDensity)); + + %% ================ output =============== + RootProperties.frac = flipud(rootFrac);% [unitless] + RootProperties.crossArea = rootCrossArea; %[m2] + RootProperties.spac = flipud(rootSpac); % [m] + RootProperties.lengthDensity = flipud(rootLengthDensity); %[m/m3] + +end \ No newline at end of file From b752bec9c6e16baef6b0f93978cac84baa38656e Mon Sep 17 00:00:00 2001 From: ZSong Date: Tue, 13 Aug 2024 17:24:32 +0200 Subject: [PATCH 05/17] remove biochemical setting --- src/STEMMUS_SCOPE.m | 2 +- src/setScenario.m | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 07af247b..e41805f6 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -72,7 +72,7 @@ RWUg = []; %% Set senarios - [biochemical, gsMethod, phwsfMethod] = setScenario(gsOption, phsOption); + [gsMethod, phwsfMethod] = setScenario(gsOption, phsOption); %% 1. define Constants Constants = io.define_constants(); diff --git a/src/setScenario.m b/src/setScenario.m index 25531713..6b831b9f 100644 --- a/src/setScenario.m +++ b/src/setScenario.m @@ -1,4 +1,4 @@ -function [biochemical, gsMethod, phwsfMethod] = setScenario(gsOption, phsOption) +function [gsMethod, phwsfMethod] = setScenario(gsOption, phsOption) %{ This function is used to set the stomatal conductance scheme and plant hydraulic pathway option. @@ -20,7 +20,7 @@ %} %% Set scenario - biochemical = @biochemical; + %biochemical = @biochemical; if strcmp(gsOption, 'BallBerry') gsMethod = 1; elseif strcmp(gsOption, 'Medlyn') From ef611c3e7a0dca6d1bae1554959787dbcb627086 Mon Sep 17 00:00:00 2001 From: ZSong Date: Tue, 13 Aug 2024 22:14:41 +0200 Subject: [PATCH 06/17] plant water potential --- src/calPlantWaterPotential.m | 198 +++++++++++++++++++++++++++++++++++ 1 file changed, 198 insertions(+) create mode 100644 src/calPlantWaterPotential.m diff --git a/src/calPlantWaterPotential.m b/src/calPlantWaterPotential.m new file mode 100644 index 00000000..933388f5 --- /dev/null +++ b/src/calPlantWaterPotential.m @@ -0,0 +1,198 @@ +function [psiLeaf, psiStem, psiRoot, kSoil2Root, kRoot2Stem, kStem2Leaf, phwsfLeaf, TempVar] = calPlantWaterPotential(Trans,Ks, Ksoil, ParaPlant,... + RootProperties, soilDepthB2T, lai, sfactor, psiSoil, canopyHeight, bbx, TestPHS, iTimeStep) +% Calculation of plant hydraulic conductance among plant components + +% Input: +% Trans: Transpiration [m/s] +% Ks : saturated soil hydraulic conductivity [m s-1] +% Ksoil: unsaturated soil hydraulic conductivity [m s-1] +% ParaPlant: A structure contains plant parameters +% RootProperties: A structure contains root properties +% soilDepthB2T: An array contains soil depth of each soil layer to soil +% surface. (direction: from bottom to top) +% +% lai: An array contains LAI +% sfactor: soil water stress factor, WSF + +% Output: +% psiLeaf: leaf water potential [m] +% psiStem: stem water potential [m] +% psiRoot: root water potential [m] +% kSoil2Root: hydraulic conductance from soil to root +% kRoot2Stem: hydraulic conductance from root to stem +% kStem2Leaf: hydraulic conductance from stem to leaf +% phwsfLeaf : leaf water stress factor + + %% +++++++++++++++++++++++++ PHS ++++++++++++++++++++++++++++++++ + % ---------------------------------------------------------------- + % qSoil2Root = kSoil2Root * (psiSoil - psiRoot - soilDepth) + % + % kSoil2Root * (psiSoil - soilDepth) - qSoil2Root + % psiRoot = -------------------------------------------------- + % kSoil2Root + + % qRoot2Stem = kRoot2Stem * SAI * (psiRoot - psiStem - canopyHeight) + % qRoot2Stem + % psiStem = psiRoot - canopyHeight - ----------------- + % kRoot2Stem * SAI + % + + % qStem2Leaf = kStem2Leaf * LAI * (psiLeaf - psiStem) + % qStem2Leaf + % psiLeaf = psiStem - --------------------- + % kStem2Leaf *LAI + % + % ---------------------------------------------------------------- + + %% =============== Input variables ================ + Krootmax = ParaPlant.Krootmax; + Kstemmax = ParaPlant.Kstemmax; + Kleafmax = ParaPlant.Kleafmax; + p50Root = ParaPlant.p50Root; + p50Stem = ParaPlant.p50Stem; + p50Leaf = ParaPlant.p50Leaf; + shapeFactorRoot = ParaPlant.shapeFactorRoot; + shapeFactorStem = ParaPlant.shapeFactorStem; + shapeFactorLeaf = ParaPlant.shapeFactorLeaf; + rootLateralLength = ParaPlant.rootLateralLength; + s2l = ParaPlant.s2l; + phwsfMethod = ParaPlant.phwsfMethod; + + rootSpac = RootProperties.spac; + rootFrac = RootProperties.frac; + lengthDensity = RootProperties.lengthDensity; + + endOfSeason = TestPHS.endOfSeason; + % inverse soilDepth +% soilDepthflip = flipud(soilDepth); + + % Q_soil2root = Q_root2stem = Q_stem2leaf = Transpiration + qSoil2Root = Trans; % unit [m/s] + qRoot2Stem = Trans; + qStem2Leaf = Trans; + + numSoilLayer = length(Ksoil); +% %% =================== reverse soil variables ============== +% KsoilR = flipud(Ksoil); +% KsR = flipud(Ks); + %% =================== root area index ===================== + % new fine-root carbon to new foliage carbon ratio + if iTimeStep <= endOfSeason % during growth season + froot2leaf=0.3*3*exp(-0.15*lai)/(exp(-0.15*lai)+2*sfactor); + if froot2leaf<0.15 + froot2leaf=0.15; + end + if froot2leaf>0.3 + froot2leaf = 0.3; + end + else % when leaves start desiccation, set carbon allocation for roots growth as a constant + froot2leaf = TestPHS.froot2leaf; + end + froot2leaf=max(0.001, froot2leaf); + + % stem area index + sai = s2l .*lai; + + % root area index + % rai = (sai + lai) * f_{root-shoot} * r_i (Kennedy et al. 2019 JAMES) + rai = (sai + lai).* rootFrac .* froot2leaf; + + %% =================== soil to root conductance ===================== + soilConductance = min(Ks' , Ksoil) ./100 ./ rootSpac ; % 100 is a transfer factor from [cm/s] to [m/s] + + phwsfRoot = PlantHydraulicsStressFactor(psiSoil, p50Root, shapeFactorRoot, phwsfMethod); + phwsfRoot = max(phwsfRoot, 1e-2); + + rootConductance = phwsfRoot .* rai .* Krootmax./(rootLateralLength + soilDepthB2T./100); % unit [m/s] + + soilConductance = max(soilConductance, 1e-16); + rootConductance = max(rootConductance, 1e-16); + kSoil2Root = 1./(1./soilConductance + 1./rootConductance); % unit [m/s] + + %% =================== root water potential ======================== + % Q_soil2root = Q_root2stem = Q_stem2leaf = Transpiration + if (abs(sum(kSoil2Root,1)) == 0) + % for saturated condition + psiRoot = sum((psiSoil - soilDepthB2T./100)) / numSoilLayer; + else + % for unsaturated condition +% psiRoot = (sum(kSoil2Root.*(psiSoil - soilDepthB2T./100).*bbx) - qSoil2Root) / sum(kSoil2Root.*bbx); + + % %-consider the root fraction on calculation of root water potential + % -------------------rootFrac------------------------------ +% rootFracFactor = (rootFrac.*bbx)./sum(rootFrac.*bbx); +% psiRoot_temp = ((kSoil2Root.*(psiSoil - soilDepthB2T./100)) - qSoil2Root.*rootFracFactor) ./ (kSoil2Root); +% psiRoot = nansum(psiRoot_temp.*rootFracFactor ,1); + % --------------------------------------------------------- + + % consider the root fraction via root length density + % ------------------ rootLengthDensityFrac ----------------- + rootLengthDensityFrac = (lengthDensity.*bbx)./sum(lengthDensity.*bbx); + psiRoot_temp = ((kSoil2Root.*(psiSoil - soilDepthB2T./100)) - qSoil2Root.*rootLengthDensityFrac) ./ (kSoil2Root); + psiRoot = nansum(psiRoot_temp.*rootLengthDensityFrac ,1); + + if ~isreal(psiRoot) +% psiRoot = sum(psiSoil.*bbx)/sum(bbx); % root zone averaged soil water potential + psiRoot = sum(psiSoil.*rootLengthDensityFrac) + end + psiRoot = max(psiRoot, -1000); % -10MPa to m = -1019.368 + psiRoot = min(0, psiRoot); + end + + %% =================== stem water potential ======================== + phwsfStem = PlantHydraulicsStressFactor(psiRoot, p50Stem, shapeFactorStem, phwsfMethod); + phwsfStem = max(phwsfStem, 0.1); + if ( sai>0 && phwsfStem >0 ) + % stem hydraulic conductance + kRoot2Stem = ParaPlant.Kstemmax ./ canopyHeight .* phwsfStem; + psiStem = psiRoot - canopyHeight - qRoot2Stem ./sai ./ kRoot2Stem; + else + psiStem = psiRoot - canopyHeight; + end + + %% ===================== leaf water potential ==================== + phwsfStem2Leaf = PlantHydraulicsStressFactor(psiStem, p50Leaf, shapeFactorLeaf, phwsfMethod); + if (lai>0 && phwsfStem2Leaf>0) + % leaf hydraulic conductance + kStem2Leaf = ParaPlant.Kleafmax .* phwsfStem2Leaf; + + % leaf water potential + psiLeaf = psiStem - qStem2Leaf ./ lai ./kStem2Leaf; + else + psiLeaf = psiStem; + end + phwsfLeaf = PlantHydraulicsStressFactor(psiLeaf, -150, shapeFactorLeaf, phwsfMethod); + + %% ==================== set complex value ====================== + if ~isreal(psiRoot) +% psiRoot = sum(psiSoil.*bbx)/sum(bbx); % root zone averaged soil water potential + psiRoot = sum(psiSoil.*rootFracFactor) + end + if ~isreal(psiStem) + psiStem = psiRoot; + end + + + %% ================== set default conductance ================== + if ~exist('kSoil2Root','var') + kSoil2Root = NaN; + end + + if ~exist('kRoot2Stem','var') + kRoot2Stem = NaN; + end + + if~exist('kStem2Leaf', 'var') + kStem2Leaf = NaN; + end + + %% ======================= Temp ================================== + TempVar.froot2leaf = froot2leaf; + TempVar.sai = sai; + TempVar.rai = rai; + TempVar.phwsfRoot = phwsfRoot; + TempVar.phwsfStem2Leaf = phwsfStem2Leaf; + TempVar.soilConductance = soilConductance; + TempVar.rootConductance = rootConductance; + +end \ No newline at end of file From 1dc0929b5dea4704225625c43c49b252a0141a65 Mon Sep 17 00:00:00 2001 From: ZSong Date: Tue, 13 Aug 2024 22:15:03 +0200 Subject: [PATCH 07/17] plant water stress factor --- src/PlantHydraulicsStressFactor.m | 96 +++++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) create mode 100644 src/PlantHydraulicsStressFactor.m diff --git a/src/PlantHydraulicsStressFactor.m b/src/PlantHydraulicsStressFactor.m new file mode 100644 index 00000000..7c8a23eb --- /dev/null +++ b/src/PlantHydraulicsStressFactor.m @@ -0,0 +1,96 @@ +% psi50 = ParaPlant.psi50; +% ck = ParaPlant.ck; +function phwsf = PlantHydraulicsStressFactor(psi, psi50, shapeFactor, phwsfMethod) +% This function is to calculate water stress factor based on plant +% hydraulisc theory. +% Input: +% psi: water potential +% psi50: parameter of plant hydraulic pathway +% shapeFactor: a struct contains parameter of plant hydraulic pathway +% ck for CLM5 +% a for ED2 +% m for PHS +% phwsf_method : phwsf_method, the default approah is weibull method +% +% Output: +% phwsf: plant hydraulic water stress factor, scale from 0 to 1. +% The value of 0 means that a severe water stress, and the value of +% 1 means there is no water stress. +% +% References: +% 1. D. Kennedy et al_2019_JAMESM_Implementing Plant Hydraulics in the Community Land Model, Version 5, DOI: https://doi.org/10.1029/2018MS001500 +% 2. X. Xu et al_2016_New Phytol_Diversity in plant hydraulic traits explains seasonal and inter-annual variations of vegetation dynamics in seasonally dry tropical forests, DOI: 10.1111/nph.14009 + + + % ========== define phwsf method =============== + if nargin < 4 + phwsfMethod = 'CLM5'; + else + phwsfMethod = phwsfMethod; + end + + % ========= calculate phwsf =================== + switch phwsfMethod + case 'CLM5' + ckCLM = shapeFactor; % shape factor in CLM5 scheme + phwsf = CLM5(psi, psi50, ckCLM); + phwsf(phwsf < 5e-5) = 0; + case 'ED2' + aED2 = shapeFactor; % shape factor in ED2 scheme + phwsf = ED2(psi, psi50, aED2); + case 'STEMMUS-SCOPE' + mPHS = shapeFactor; % shape factor in STEMMUS-SCOPE-PHS + p0 = -0.33; + phwsf = PHS(psi, psi50, mPHS) +% case '' + otherwise + phwsf = NaN; + fprintf('phwsf method need to be defined\n.') + end +end + +function phwsf = CLM5(psi, psi50, ck) +%{ + This function calculated plant hydraulic water stress factor based on + the scheme of CLM5 + Input: + psi : leaf water potential (m) + psi50 : P50 value (m) + ck : shape factor (-) + Output: + phwsf : plant hydraulic water stress factor, range from 0 to 1. +%} + phwsf = 2.^(-(psi./psi50)).^ck; + phwsf(phwsf < 5e-5) = 0; + +end + +function phwsf = ED2(psi, psi50, a) +%{ + This function calculated plant water stress factor based on the scheme + of ED2. (Xu, Xiangtao_2016_new phytologist) +%} + + phwsf = (1+(psi./psi50).^a).^-1; +end + + +function phwsf = PHS(psi, psi50, m) +%{ + A new plant water stress function format based on soil water stress + factor. Since the observed plant water potential with the unit of MPa, + thus, we translate the unit from m to MPa at first. Then use curve + fitting to retrive the parameters. +%} + m2MPa = 1*9810/1e6; + psiMPa = psi .* m2MPa; + psi50MPa = psi50 .* m2MPa; + psi0 = -0.33; % psi0 is p0, we use the soil water potential at the field capacity to represent this value. + % -m .* psi0 = psi0 + (psi1.5MPa + psi0)/2 + phwsf = (1+exp(-m .* psi0 * (psiMPa - psi50MPa))).^-1; +% % -m .* psi0 = psi0 + (psi1.5MPa + psi0)/2 +end + + + + From b644df55941ce31bc1489c81cd69d6758cffbe6a Mon Sep 17 00:00:00 2001 From: ZSong Date: Tue, 13 Aug 2024 22:16:22 +0200 Subject: [PATCH 08/17] plant hydraulics --- src/STEMMUS_SCOPE.m | 10 +- src/biochemical.m | 1 + src/calc_rsoil.m | 10 +- src/ebal.m | 229 +++++++++++++++++++++++++++++++++++--------- src/heatfluxes.m | 4 +- 5 files changed, 204 insertions(+), 50 deletions(-) diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index e41805f6..b33679a0 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -89,6 +89,7 @@ [RootProperties,soilDepth] = calRootProperties(SiteProperties, ParaPlant, numSoilLayer, soilThickness, RTB); Rl = RootProperties.lengthDensity; Ztot = soilDepth; + soilDepthB2T = flipud(soilDepth); %% 2. simulation options @@ -127,6 +128,7 @@ useXLSX = 1; % set it to 1 or 0, the current stemmus-scope does not support 0 [ScopeParameters, options] = parameters.loadParameters(options, useXLSX, X, F, N); options.gsMethod = gsMethod; % 1 for BallBerry's method; 2 for Medlyn's method + options.plantHydraulics = phsOption; % Define the location information ScopeParameters.LAT = SiteProperties.latitude; % latitude @@ -332,7 +334,8 @@ % (this is to not repeat the save-workspace code). disp('Finished model initialization'); end - +TestPHS.psiLeafIni = 0-SiteProperties.canopy_height; +TestPHS.endOfSeason = find(ScopeParameters.LAI == max(ScopeParameters.LAI)); % If the runMode is update, retrieve the (possibly) updated state. % The state can be modified by the STEMMUS_SCOPE BMI. See PyStemmusScope. if strcmp(bmiMode, 'update') @@ -490,10 +493,11 @@ switch options.calc_ebal case 1 - [iter, fluxes, rad, thermal, profiles, soil, RWU, frac, WaterStressFactor, WaterPotential] ... + [iter, fluxes, rad, thermal, profiles, soil, RWU, frac, WaterStressFactor, WaterPotential, TestPHS] ... = ebal(iter, options, spectral, rad, gap, ... leafopt, angles, meteo, soil, canopy, leafbio, xyt, k, profiles, Delt_t, ... - Rl, SoilVariables, VanGenuchten, InitialValues, GroundwaterSettings); + Rl, SoilVariables, VanGenuchten, InitialValues, GroundwaterSettings,... + SiteProperties, ParaPlant, RootProperties, soilDepthB2T, TestPHS, KT); if options.calc_fluor if options.calc_vert_profiles [rad, profiles] = RTMf(spectral, rad, soil, leafopt, canopy, gap, angles, profiles); diff --git a/src/biochemical.m b/src/biochemical.m index a5f79463..bb261b3a 100644 --- a/src/biochemical.m +++ b/src/biochemical.m @@ -553,6 +553,7 @@ biochem_out.Fm_Fo = fm ./ fo; % parameters used for curve fitting biochem_out.Ft_Fo = fs ./ fo; % parameters used for curve fitting biochem_out.qQ = qQ; + biochem_out.VPD_l2b = mean(VPD_l2b,'all'); return end % end of function biochemical diff --git a/src/calc_rsoil.m b/src/calc_rsoil.m index a37cc034..89b437e2 100644 --- a/src/calc_rsoil.m +++ b/src/calc_rsoil.m @@ -1,13 +1,19 @@ -function [PSIs, rsss, rrr, rxx] = calc_rsoil(Rl, ModelSettings, SoilVariables, VanGenuchten, bbx, GroundwaterSettings) +function [PSIs, rsss, rrr, rxx, psiSoilTemp, Ksoil] = calc_rsoil(Rl, ModelSettings, SoilVariables, VanGenuchten, bbx, GroundwaterSettings) % load Constants Constants = io.define_constants(); SMC = SoilVariables.Theta_LL(1:end - 1, 2); Se = (SMC - VanGenuchten.Theta_r') ./ (VanGenuchten.Theta_s' - VanGenuchten.Theta_r'); + if any(Se<0) + Se(Se<0)=0; + end + Ksoil = SoilVariables.Ks' .* Se.^Constants.l .* (1 - (1 - Se.^(1 ./ VanGenuchten.m')).^(VanGenuchten.m')).^2; if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling - PSIs = -((Se.^(-1 ./ VanGenuchten.m') - 1).^(1 ./ VanGenuchten.n')) ./ (VanGenuchten.Alpha * 100)' .* bbx; + psiSoilTemp = -((Se.^(-1 ./ VanGenuchten.m') - 1).^(1 ./ VanGenuchten.n')) ./ (VanGenuchten.Alpha * 100)'; + psiSoilTemp = max(psiSoilTemp, -5e3); + PSIs = psiSoilTemp.* bbx; else % Groundwater coupling is activated, added by Mostafa (see https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/231) % Change PSIs with SoilVariables.hh to correct hydraulic head (matric potential + gravity) of the saturated layers for i = 1:ModelSettings.NL diff --git a/src/ebal.m b/src/ebal.m index e29eff79..bb91af94 100644 --- a/src/ebal.m +++ b/src/ebal.m @@ -1,7 +1,8 @@ -function [iter, fluxes, rad, thermal, profiles, soil, RWU, frac, WaterStressFactor, WaterPotential] ... +function [iter, fluxes, rad, thermal, profiles, soil, RWU, frac, WaterStressFactor, WaterPotential, TestPHS] ... = ebal(iter, options, spectral, rad, gap, leafopt, ... angles, meteo, soil, canopy, leafbio, xyt, k, profiles, Delt_t, ... - Rl, SoilVariables, VanGenuchten, InitialValues, GroundwaterSettings) + Rl, SoilVariables, VanGenuchten, InitialValues, GroundwaterSettings, ... + SiteProperties, ParaPlant, RootProperties, soilDepthB2T, TestPHS, KT) %{ function ebal.m calculates the energy balance of a vegetated surface @@ -109,6 +110,7 @@ resistance of leaves (or biochemical_MD12: alternative) Ca = meteo.Ca; Ts = soil.Ts; p = meteo.p; + % RH = meteo.RH; if options.soil_heat_method < 2 && options.simulation == 1 if k > 1 Deltat = Delt_t; % Duration of the time interval (s) @@ -176,9 +178,31 @@ resistance of leaves (or biochemical_MD12: alternative) PSI = 0; [bbx] = Max_Rootdepth(InitialValues.bbx); - [PSIs, rsss, rrr, rxx] = calc_rsoil(Rl, ModelSettings, SoilVariables, VanGenuchten, bbx, GroundwaterSettings); + [psiSoil, rsss, rrr, rxx, TestPHS.psiSoilAll(:, KT), Ksoil] = calc_rsoil(Rl, ModelSettings, SoilVariables, VanGenuchten, bbx, GroundwaterSettings); + TestPHS.rsssTot(:, KT) = rsss; + TestPHS.rrrTot(:, KT) = rrr; + TestPHS.rxxTot(:, KT) = rxx; + [sfactor] = calc_sfactor(Rl, VanGenuchten.Theta_s, VanGenuchten.Theta_r, SoilVariables.Theta_LL, bbx, Ta, VanGenuchten.Theta_f); - PSIss = PSIs(ModelSettings.NL, 1); + + PSIss = psiSoil(ModelSettings.NL, 1); + + % initial leaf water potental = soil water potential - gravitational potential + canopyHeight = SiteProperties.canopy_height(KT); + + % psiLeaf = 0-canopyHeight; + psiLeaf = TestPHS.psiLeafIni(KT); + PSI = 0; + % psiAir = air_water_potential(RH, Ta); + airPress_m = meteo.p .*1e2 ./9810; + airPress_hPa = meteo.p; + + % options.plantHydraulics = 1; % Indicating whether to use PHS: 1 PHS open; 0 PHS close. + if options.plantHydraulics + phwsf = PlantHydraulicsStressFactor(psiLeaf, ParaPlant.p50Leaf, ParaPlant.shapeFactorLeaf, ParaPlant.phwsfMethod); + else + phwsf = 1; + end %% 2. Energy balance iteration loop % 'Energy balance loop (Energy balance and radiative transfer) @@ -281,6 +305,8 @@ resistance of leaves (or biochemical_MD12: alternative) rcwh = biochem_out.rcw; qEh = biochem_out.qE; % vCaemmerer- Magnani does not generate this parameter (dummy value) Knh = biochem_out.Kn; + VPDh = biochem_out.VPD_l2b; + % for sunlit leaves biochem_in.T = Tcu; @@ -304,6 +330,7 @@ resistance of leaves (or biochemical_MD12: alternative) rcwu = biochem_out.rcw; qEu = biochem_out.qE; Knu = biochem_out.Kn; + VPDu = biochem_out.VPD_l2b; Pinh = rad.Pnh; Pinu = rad.Pnu; @@ -315,51 +342,146 @@ resistance of leaves (or biochemical_MD12: alternative) % 2.4. Fluxes (latent heat flux (lE), sensible heat flux (H) and soil heat flux G % in analogy to Ohm's law, for canopy (c) and soil (s). All in units of [W m-2] - % soil.PSIs; + % soil.psiSoil; rss = soil.rss; rac = (LAI + 1) * (raa + rawc); ras = (LAI + 1) * (raa + raws); - for i = 1:30 - [lEch, Hch, ech, Cch, lambdah, sh] = heatfluxes(rac, rcwh, Tch, ea, Ta, e_to_q, PSI, Ca, Cih, es_fun, s_fun); - [lEcu, Hcu, ecu, Ccu, lambdau, su] = heatfluxes(rac, rcwu, Tcu, ea, Ta, e_to_q, PSI, Ca, Ciu, es_fun, s_fun); - [lEs, Hs, ~, ~, lambdas, ss] = heatfluxes(ras, rss, Ts, ea, Ta, e_to_q, PSIss, Ca, Ca, es_fun, s_fun); - - % if any( ~isreal( Cch )) || any( ~isreal( Ccu(:) )) - % error('Heatfluxes produced complex values for CO2 concentration!') - % end - - % if any( Cch < 0 ) || any( Ccu(:) < 0 ) - % error('Heatfluxes produced negative values for CO2 concentration!') - % end - - % integration over the layers and sunlit and shaded fractions - Hstot = Fs * Hs; - Hctot = LAI * (Fc * Hch + equations.meanleaf(canopy, Hcu, 'angles_and_layers', Ps)); - Htot = Hstot + Hctot; - %%%%%% Leaf water potential calculate - lambda1 = (2.501 - 0.002361 * Ta) * 1E6; - lEctot = LAI * (Fc * lEch + equations.meanleaf(canopy, lEcu, 'angles_and_layers', Ps)); % latent heat leaves - if isreal(lEctot) && lEctot < 1000 && lEctot > -300 - else - lEctot = 0; - end - Trans = lEctot / lambda1 / 1000; % unit: m s-1 - AA1 = PSIs ./ (rsss + rrr + rxx); - AA2 = 1 ./ (rsss + rrr + rxx); - BB1 = AA1(~isnan(AA1)); - BB2 = AA2(~isinf(AA2)); - PSI1 = (sum(BB1) - Trans) / sum(BB2); - if isnan(PSI1) - PSI1 = -1; + + if options.plantHydraulics + % ======================== PHS open =========================== + for i=1:30 + [lEch,Hch,ech,eih, Cch,lambdah,sh, delta_eh, delta_th] = heatfluxes(rac,rcwh,Tch,ea,Ta,e_to_q,psiLeaf,Ca,Cih,es_fun,s_fun); + [lEcu,Hcu,ecu,eiu, Ccu,lambdau,su, delta_eu, delta_tu] = heatfluxes(rac,rcwu,Tcu,ea,Ta,e_to_q,psiLeaf,Ca,Ciu,es_fun,s_fun); + [lEs,Hs,~,~,~,lambdas,ss, delta_es, delta_ts] = heatfluxes(ras,rss,Ts ,ea,Ta,e_to_q,PSIss,Ca,Ca,es_fun,s_fun); + + + % integration over the layers and sunlit and shaded fractions + Hstot = Fs*Hs; + Hctot = LAI*(Fc*Hch + equations.meanleaf(canopy,Hcu,'angles_and_layers',Ps)); + Htot = Hstot + Hctot; + + %%%%%% Leaf water potential calculate + lambda1 = (2.501-0.002361*Ta)*1E6; + lEctot = LAI*(Fc*lEch + equations.meanleaf(canopy,lEcu,'angles_and_layers',Ps)); % latent heat leaves + if (isreal(lEctot) && lEctot<1000 && lEctot>-300) + else + lEctot=0; + end + Trans = lEctot/lambda1/1000; % total canopy transpiration: unit: m s-1 + % Trans_t = lEct .* LAI./lambda1./1000; + + [psiLeaf_temp, psiStem, psiRoot, kSoil2Root, kRoot2Stem, kStem2Leaf, phwsf, TempVar] = calPlantWaterPotential(Trans,SoilVariables.Ks, ... + Ksoil, ParaPlant, RootProperties, soilDepthB2T, LAI, sfactor, psiSoil, canopyHeight, bbx, TestPHS, KT); + + if isnan(psiLeaf_temp)|~isreal(psiLeaf_temp) + psiLeaf_temp = -1; + end + + if abs(psiLeaf - psiLeaf_temp)<0.01 + break + end + psiLeaf = 0.5 * (psiLeaf + psiLeaf_temp); end - if ~isreal(PSI1) - PSI1 = -1; + + % if phwsf is a complex value, set it as sfactor + if ~isreal(phwsf) + phwsf = sfactor; end - if abs(PSI - PSI1) < 0.01 - break + + % stomatal resistance + canopyStoResis = Fc*rcwh + equations.meanleaf(canopy,rcwu,'angles_and_layers',Ps); + + % canopy conductance = 1/(stomatal resistance + aerodynamic resistance) + canopyConduct = 1./(canopyStoResis + rac); + + phsTrans = canopyConduct .* LAI .* VPDu./airPress_hPa; + %% ====================== root water uptake ===================== + rootWaterUptake = kSoil2Root .* (psiSoil - psiRoot - soilDepthB2T./100).*bbx; + + + TestPHS.psiStemTot(KT) = psiStem; + TestPHS.psiRootTot(KT) = psiRoot; + TestPHS.psiSoilTot(:,KT) = psiSoil; % psiSoil + TestPHS.psiSoilTotMean(KT) = sum(psiSoil.*bbx)/sum(bbx); + TestPHS.psiLeafTot(KT) = psiLeaf; + TestPHS.kSoil2RootTot(:,KT) = kSoil2Root; + TestPHS.kSoil2RootTotMean(KT) = sum(kSoil2Root .* bbx)/sum(bbx); + TestPHS.kRoot2StemTot(KT) = kRoot2Stem; + TestPHS.kStem2LeafTot(KT) = kStem2Leaf; + TestPHS.phwsfTot(KT) = phwsf; + TestPHS.transTot(KT) = Trans; + + % TestPHS.psiAirTot(KT) = psiAir; + % TestPHS.kLeaf2AirTot(KT) = Trans./(psiLeaf - psiAir); + + % TestPHS.trans_t(KT) = Trans_t; + TestPHS.phsTrans(KT) = phsTrans; + TestPHS.canopyStoResisTot(KT) = canopyStoResis; + TestPHS.racTot(KT) = rac; + TestPHS.LAI(KT) = LAI; + TestPHS.canopyConductTot(KT) = canopyConduct; + + TestPHS.froot2leaf = TempVar.froot2leaf; + TestPHS.froot2leafTot(KT) = TempVar.froot2leaf; + TestPHS.saiTot(KT) = TempVar.sai; + TestPHS.raiTot(:,KT) = TempVar.rai; + + TestPHS.soilConductanceTot(:,KT) = TempVar.soilConductance; + TestPHS.rootConductanceTot(:,KT) = TempVar.rootConductance; + TestPHS.phwsfRootTot(:,KT) = TempVar.phwsfRoot; + TestPHS.phwsfStem2LeafTot(KT) = TempVar.phwsfStem2Leaf; + + TestPHS.raiTotMean(KT) = sum(TempVar.rai .* bbx)/sum(bbx); + TestPHS.soilConductanceTotMean(KT) = sum(TempVar.soilConductance .* bbx)/sum(bbx); + TestPHS.rootConductanceTotMean(KT) = sum(TempVar.rootConductance .* bbx)/sum(bbx); + % ===================== PHS close ============================== + else + for i = 1:30 + [lEch, Hch, ech, ~, Cch, lambdah, sh, ~, ~] = heatfluxes(rac, rcwh, Tch, ea, Ta, e_to_q, PSI, Ca, Cih, es_fun, s_fun); + [lEcu, Hcu, ecu, ~, Ccu, lambdau, su, ~, ~] = heatfluxes(rac, rcwu, Tcu, ea, Ta, e_to_q, PSI, Ca, Ciu, es_fun, s_fun); + [lEs, Hs, ~, ~, ~, lambdas, ss,~,~] = heatfluxes(ras, rss, Ts, ea, Ta, e_to_q, PSIss, Ca, Ca, es_fun, s_fun); + + % if any( ~isreal( Cch )) || any( ~isreal( Ccu(:) )) + % error('Heatfluxes produced complex values for CO2 concentration!') + % end + + % if any( Cch < 0 ) || any( Ccu(:) < 0 ) + % error('Heatfluxes produced negative values for CO2 concentration!') + % end + + % integration over the layers and sunlit and shaded fractions + Hstot = Fs * Hs; + Hctot = LAI * (Fc * Hch + equations.meanleaf(canopy, Hcu, 'angles_and_layers', Ps)); + Htot = Hstot + Hctot; + %%%%%% Leaf water potential calculate + lambda1 = (2.501 - 0.002361 * Ta) * 1E6; + lEctot = LAI * (Fc * lEch + equations.meanleaf(canopy, lEcu, 'angles_and_layers', Ps)); % latent heat leaves + if isreal(lEctot) && lEctot < 1000 && lEctot > -300 + else + lEctot = 0; + end + Trans = lEctot / lambda1 / 1000; % unit: m s-1 + AA1 = psiSoil ./ (rsss + rrr + rxx); + AA2 = 1 ./ (rsss + rrr + rxx); + BB1 = AA1(~isnan(AA1)); + BB2 = AA2(~isinf(AA2)); + psiLeaf_temp = (sum(BB1) - Trans) / sum(BB2); + if isnan(psiLeaf_temp) | ~isreal(psiLeaf_temp) + psiLeaf_temp = -1; + end + % if ~isreal(PSI1) + % PSI1 = -1; + % end + if abs(psiLeaf - psiLeaf_temp) < 0.01 + break + end + psiLeaf = 0.5 .* (psiLeaf + psiLeaf_temp); end - PSI = (PSI + PSI1) / 2; + TestPHS.psiSoilTot(:, KT) = psiSoil; + TestPHS.psiSoilTotMean(KT) = sum(psiSoil.*bbx)/sum(bbx); + TestPHS.psiLeafTot(KT) = psiLeaf; end + PSItot(KT) = psiLeaf; %%%%%%% if SoilHeatMethod == 2 @@ -436,6 +558,19 @@ resistance of leaves (or biochemical_MD12: alternative) Lot_ = equations.Planck(spectral.wlS', Tbr); rad.LotBB_ = Lot_; % Note that this is the blackbody radiance! + %% =============== debug: output resistance 20221205 Z.So ========== + TestPHS.rssTot(KT) = rss; % Surface resistance of soil for vapour transport + TestPHS.racTot(KT) = rac; % aerodynamic resistance for heat in canopy + TestPHS.rasTot(KT) = ras; % aerodynamic resistance for heat in soil + TestPHS.rcwhTot(:,KT) = rcwh; % stomatal resistance of sunlit leaf + TestPHS.rcwuTot(:,KT) = reshape(rcwu,[],1); % stomatal resistance of sunlit leaf + TestPHS.rcwTot(KT) = Fc*rcwh+equations.meanleaf(canopy,rcwu,'angles_and_layers',Ps); % stomatal resistance of sunlit leaf + TestPHS.gamTot(KT) = GAM; + TestPHS.delta_ecTot(KT) = Fc*delta_eh+equations.meanleaf(canopy, delta_eu, 'angles_and_layers',Ps); + TestPHS.delta_tcTot(KT) = Fc*delta_th+equations.meanleaf(canopy, delta_tu, 'angles_and_layers',Ps); + TestPHS.delta_esTot(KT) = Fs*delta_es; + TestPHS.delta_tsTot(KT) = Fs*delta_ts; + %% 3. Print warnings whenever the energy balance could not be solved if counter >= maxit warning(['\n warning: maximum number of iteratations exceeded', ... @@ -552,7 +687,12 @@ resistance of leaves (or biochemical_MD12: alternative) fluxes.Au = Au; fluxes.Ah = Ah; - RWU = (PSIs - PSI) ./ (rsss + rrr + rxx) .* bbx; + if options.plantHydraulics + RWU = rootWaterUptake; + else + RWU =(psiSoil - psiLeaf)./(rsss+rrr+rxx).*bbx; + end + nn = numel(RWU); for i = 1:nn if isnan(RWU(i)) @@ -565,7 +705,7 @@ resistance of leaves (or biochemical_MD12: alternative) end end frac = RWU ./ abs(sum(sum(RWU))); - RWU = (PSIs - PSI) ./ (rsss + rrr + rxx) .* bbx; + % RWU = (PSIs - PSI) ./ (rsss + rrr + rxx) .* bbx; RWU = real(RWU); for i = 1:nn if isnan(RWU(i)) @@ -580,5 +720,6 @@ resistance of leaves (or biochemical_MD12: alternative) % return WaterStressFactor.soil = sfactor; + WaterStressFactor.plant = phwsf; WaterPotential.leaf = PSI; end diff --git a/src/heatfluxes.m b/src/heatfluxes.m index 7729f14e..2b4c33e9 100644 --- a/src/heatfluxes.m +++ b/src/heatfluxes.m @@ -1,4 +1,4 @@ -function [lE, H, ec, Cc, lambda, s] = heatfluxes(ra, rs, Tc, ea, Ta, e_to_q, PSI, Ca, Ci, es_fun, s_fun) +function [lE, H, ec, ei, Cc, lambda, s, delta_e, delta_t] = heatfluxes(ra, rs, Tc, ea, Ta, e_to_q, PSI, Ca, Ci, es_fun, s_fun) % load Constants Constants = io.define_constants(); @@ -47,3 +47,5 @@ H = (rhoa * cp) ./ ra .* (Tc - Ta); % [W m-2] Sensible heat flux ec = ea + (ei - ea) * ra ./ (ra + rs); % [W m-2] vapour pressure at the leaf surface Cc = Ca - (Ca - Ci) .* ra ./ (ra + rs); % [umol m-2 s-1] CO2 concentration at the leaf surface + delta_e = ei-ea; + delta_t = Tc-Ta; \ No newline at end of file From 9ddaad54757177a5fab32f4c8df16fa5e33b190f Mon Sep 17 00:00:00 2001 From: ZSong Date: Wed, 14 Aug 2024 12:43:53 +0200 Subject: [PATCH 09/17] set NBChB according to site propoerties --- src/+init/setBoundaryCondition.m | 5 ++++- src/STEMMUS_SCOPE.m | 2 +- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/+init/setBoundaryCondition.m b/src/+init/setBoundaryCondition.m index 5e83869d..3007805d 100644 --- a/src/+init/setBoundaryCondition.m +++ b/src/+init/setBoundaryCondition.m @@ -1,4 +1,4 @@ -function BoundaryCondition = setBoundaryCondition(SoilVariables, ForcingData, initialLandcoverClass) +function BoundaryCondition = setBoundaryCondition(SoilVariables, ForcingData, initialLandcoverClass, siteName) % Variables information for boundary condition settings % NBCh Indicator for type of surface boundary condition on mass euqation to be applied; % "1"--Specified matric head; @@ -69,6 +69,9 @@ else NBChB = 3; end + if strcmp(siteName, "CH-HTC") + NBChB = 1; + end BChB = -9e-10; if ModelSettings.Thmrlefc == 1 % NBCT: Energy Surface B.C.: diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index b33679a0..3f1c1aa8 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -276,7 +276,7 @@ SoilConstants = io.getSoilConstants(); %% The boundary condition information settings - BoundaryCondition = init.setBoundaryCondition(SoilVariables, ForcingData, SiteProperties.landcoverClass(1)); + BoundaryCondition = init.setBoundaryCondition(SoilVariables, ForcingData, SiteProperties.landcoverClass(1), SiteProperties.sitename); DSTOR = BoundaryCondition.DSTOR; DSTOR0 = BoundaryCondition.DSTOR0; RS = BoundaryCondition.RS; From 4ee3c86dca17a6fc1bacd6f4902324986fb4ecf0 Mon Sep 17 00:00:00 2001 From: ZSong Date: Wed, 14 Aug 2024 12:45:04 +0200 Subject: [PATCH 10/17] load Initial soil depth setting from external files --- src/+io/loadSoilInitialValues.m | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/+io/loadSoilInitialValues.m b/src/+io/loadSoilInitialValues.m index ebf04fdb..aabff958 100644 --- a/src/+io/loadSoilInitialValues.m +++ b/src/+io/loadSoilInitialValues.m @@ -35,12 +35,18 @@ BtmX = BtmX; Tss = Tss; - InitND1 = 5; % Unit of it is cm. These variables are used to indicated the depth corresponding to the measurement. - InitND2 = 15; - InitND3 = 60; - InitND4 = 100; - InitND5 = 200; - InitND6 = 300; + % This input of soil layer thickness and initial depth will be revised + % later via reading csv file. I read the initial soil depth variable + % from .mat file. + % ZSong (August 2024) + if ~exist("InitND1", "var") + InitND1 = 5; % Unit of it is cm. These variables are used to indicated the depth corresponding to the measurement. + InitND2 = 15; + InitND3 = 60; + InitND4 = 100; + InitND5 = 200; + InitND6 = 300; + end if SWCC == 0 InitT0 = -1.762; % -1.75estimated soil surface temperature-1.762 InitT1 = -0.662; From 59586c6b12fc1abcb9c2da5bbcd5327e026c1c70 Mon Sep 17 00:00:00 2001 From: ZSong Date: Thu, 15 Aug 2024 12:39:23 +0200 Subject: [PATCH 11/17] PHS --- src/biochemical.m | 13 +++++++++++-- src/ebal.m | 7 ++++--- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/biochemical.m b/src/biochemical.m index bb261b3a..4dc5988a 100644 --- a/src/biochemical.m +++ b/src/biochemical.m @@ -137,6 +137,8 @@ O = biochem_in.O; p = biochem_in.p; ei = biochem_in.ei; + phwsf = biochem_in.phwsf; + plantHydraulics = biochem_in.plantHydraulics; gsMethod = biochem_in.gsMethod; g1Med = biochem_in.g1Med; @@ -283,9 +285,16 @@ % jak 2014-12-04: Add TL for C3 as well, works much better with our cotton temperature dataset (A-T) if strcmpi(Type, 'C3') && ~useTLforC3 - Vcmax = Vcmax25 .* exp(log(QTVc) .* qt) ./ TH * sfactor; + Vcmax = Vcmax25 .* exp(log(QTVc) .* qt) ./ TH;% * sfactor; else - Vcmax = Vcmax25 .* exp(log(QTVc) .* qt) ./ (TL .* TH) * sfactor; + Vcmax = Vcmax25 .* exp(log(QTVc) .* qt) ./ (TL .* TH);% * sfactor; + end + + % check it PHS open + if plantHydraulics % PHS open, use phwsf + Vcmax = Vcmax .* phwsf; + else % PHS close, use sfactor + Vcmax = Vcmax .* sfactor; end % specificity (tau in Collatz e.a. 1991) diff --git a/src/ebal.m b/src/ebal.m index bb91af94..d2f3b105 100644 --- a/src/ebal.m +++ b/src/ebal.m @@ -264,11 +264,12 @@ resistance of leaves (or biochemical_MD12: alternative) biochem_in.g1Med = leafbio.g1Med; biochem_in.g0Med = leafbio.g0Med; biochem_in.gsMethod = options.gsMethod; + biochem_in.O = meteo.Oa; biochem_in.Rdparam = leafbio.Rdparam; - biochem_in.g1Med = leafbio.g1Med; - biochem_in.g0Med = leafbio.g0Med; - + biochem_in.phwsf = phwsf; + biochem_in.plantHydraulics = options.plantHydraulics; + if options.Fluorescence_model == 2 % specific for the v.Caemmerer-Magnani model biochem_in.Tyear = leafbio.Tyear; From 3ebe43880d9bcdab1c88d455adb00111e39202aa Mon Sep 17 00:00:00 2001 From: ZSong Date: Thu, 15 Aug 2024 12:40:27 +0200 Subject: [PATCH 12/17] change default setting of BtmT --- src/+io/loadSoilInitialValues.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/+io/loadSoilInitialValues.m b/src/+io/loadSoilInitialValues.m index aabff958..38b57e8e 100644 --- a/src/+io/loadSoilInitialValues.m +++ b/src/+io/loadSoilInitialValues.m @@ -78,7 +78,8 @@ if nanmean(Ta_msr) < 0 BtmT = 0; % 9 8.1 else - BtmT = nanmean(Ta_msr); + BtmT = 14.5; + % BtmT = nanmean(Ta_msr); end if InitX0 > SaturatedMC(1) || InitX1 > SaturatedMC(1) || InitX2 > SaturatedMC(2) || ... InitX3 > SaturatedMC(3) || InitX4 > SaturatedMC(4) || InitX5 > SaturatedMC(5) || InitX6 > SaturatedMC(6) From b2b4ab7ef3dcafb555805bc705ff4db4b6cf497d Mon Sep 17 00:00:00 2001 From: Crystal-szj Date: Fri, 23 Aug 2024 10:32:24 +0200 Subject: [PATCH 13/17] revise soil parameters initialization via IS --- src/+init/soilHeteroSubroutine.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/+init/soilHeteroSubroutine.m b/src/+init/soilHeteroSubroutine.m index f9d45937..c44a833c 100644 --- a/src/+init/soilHeteroSubroutine.m +++ b/src/+init/soilHeteroSubroutine.m @@ -50,12 +50,12 @@ j = i - 1; end - SoilVariables.IS(i) = indexOfSoilType; + SoilVariables.IS(j) = indexOfSoilType; if subRoutine == 4 SoilVariables.IS(5:8) = 5; end - J = SoilVariables.IS(i); + J = SoilVariables.IS(j); [SoilVariables, VanGenuchten] = init.updateSoilVariables(SoilVariables, VanGenuchten, SoilProperties, j, J); SoilVariables.Imped(i) = ImpedF(J); From b644f6ecd546bc85c6c88701d51948c0386b7411 Mon Sep 17 00:00:00 2001 From: Crystal-szj Date: Fri, 23 Aug 2024 10:32:49 +0200 Subject: [PATCH 14/17] remove ei from output of heatflux.m --- src/ebal.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ebal.m b/src/ebal.m index d2f3b105..7596dcf2 100644 --- a/src/ebal.m +++ b/src/ebal.m @@ -351,9 +351,9 @@ resistance of leaves (or biochemical_MD12: alternative) if options.plantHydraulics % ======================== PHS open =========================== for i=1:30 - [lEch,Hch,ech,eih, Cch,lambdah,sh, delta_eh, delta_th] = heatfluxes(rac,rcwh,Tch,ea,Ta,e_to_q,psiLeaf,Ca,Cih,es_fun,s_fun); - [lEcu,Hcu,ecu,eiu, Ccu,lambdau,su, delta_eu, delta_tu] = heatfluxes(rac,rcwu,Tcu,ea,Ta,e_to_q,psiLeaf,Ca,Ciu,es_fun,s_fun); - [lEs,Hs,~,~,~,lambdas,ss, delta_es, delta_ts] = heatfluxes(ras,rss,Ts ,ea,Ta,e_to_q,PSIss,Ca,Ca,es_fun,s_fun); + [lEch,Hch,ech, Cch,lambdah,sh, delta_eh, delta_th] = heatfluxes(rac,rcwh,Tch,ea,Ta,e_to_q,psiLeaf,Ca,Cih,es_fun,s_fun); + [lEcu,Hcu,ecu, Ccu,lambdau,su, delta_eu, delta_tu] = heatfluxes(rac,rcwu,Tcu,ea,Ta,e_to_q,psiLeaf,Ca,Ciu,es_fun,s_fun); + [lEs,Hs,~,~,lambdas,ss, delta_es, delta_ts] = heatfluxes(ras,rss,Ts ,ea,Ta,e_to_q,PSIss,Ca,Ca,es_fun,s_fun); % integration over the layers and sunlit and shaded fractions From 0eadb433ef80cf2cfca4ade5780e62c97004b278 Mon Sep 17 00:00:00 2001 From: Crystal-szj Date: Fri, 23 Aug 2024 12:54:48 +0200 Subject: [PATCH 15/17] using the same parameters to calculate GAM --- src/+equations/Soil_Inertia1.m | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/+equations/Soil_Inertia1.m b/src/+equations/Soil_Inertia1.m index ce3499f7..ff08c3e0 100644 --- a/src/+equations/Soil_Inertia1.m +++ b/src/+equations/Soil_Inertia1.m @@ -7,7 +7,7 @@ Sr = SMC / theta_s; % fss = 0.58; %(sand fraction) - gamma_s = 0.27; % (soil texture dependent parameter) + gamma_s = 0.96; % (soil texture dependent parameter) dels = 1.33; % (shape parameter) ke = exp(gamma_s * (1 - power(Sr, gamma_s - dels))); @@ -15,10 +15,12 @@ phis = theta_s0; % (phis == theta_s) lambda_d = -0.56 * phis + 0.51; - QC = 0.20; % (quartz content) + QC = 0.60; % (quartz content) lambda_qc = 7.7; % (thermal conductivity of quartz, constant) - lambda_s = (lambda_qc^(QC)) * lambda_d^(1 - QC); + lambda_o = 2.0; + + lambda_s = (lambda_qc^(QC)) * lambda_o^(1 - QC); lambda_wtr = 0.57; % (thermal conductivity of water, W/m.K, constant) lambda_w = (lambda_s^(1 - phis)) * lambda_wtr^(phis); From b6c06b9b16b05608cd98504d602cfaea86ddbb20 Mon Sep 17 00:00:00 2001 From: Crystal-szj Date: Fri, 23 Aug 2024 12:55:48 +0200 Subject: [PATCH 16/17] rm ei from the output list of heatflux --- src/heatfluxes.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/heatfluxes.m b/src/heatfluxes.m index 2b4c33e9..5ae7511c 100644 --- a/src/heatfluxes.m +++ b/src/heatfluxes.m @@ -1,4 +1,4 @@ -function [lE, H, ec, ei, Cc, lambda, s, delta_e, delta_t] = heatfluxes(ra, rs, Tc, ea, Ta, e_to_q, PSI, Ca, Ci, es_fun, s_fun) +function [lE, H, ec, Cc, lambda, s, delta_e, delta_t] = heatfluxes(ra, rs, Tc, ea, Ta, e_to_q, PSI, Ca, Ci, es_fun, s_fun) % load Constants Constants = io.define_constants(); From c477a4041816ae0764b51e833ff7924d87caa6ea Mon Sep 17 00:00:00 2001 From: Crystal-szj Date: Fri, 4 Oct 2024 21:56:31 +0200 Subject: [PATCH 17/17] set gsMethod and phsOption via input_data.xlsx --- src/STEMMUS_SCOPE.m | 26 ++++++++++++++++---------- src/setScenario.m | 25 ++++++++++++++++++------- 2 files changed, 34 insertions(+), 17 deletions(-) diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 74ee08be..6e53b6ab 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -40,9 +40,6 @@ % Read the configPath file. Due to using MATLAB compiler, we cannot use run(CFG) disp (['Reading config from ', CFG]); [InputPath, OutputPath, InitialConditionPath] = io.read_config(CFG); - gsOption = 'Medlyn'; % 1 for BallBerry's method; 2 for Medlyn's method - phsOption = 2; % 1 for CLM5, 2 for ED2, 3 for PHS - % Prepare forcing and soil data [SiteProperties, SoilProperties, TimeProperties] = io.prepareInputData(InputPath); @@ -88,12 +85,11 @@ RWUs = []; RWUg = []; - %% Set senarios - [gsMethod, phwsfMethod] = setScenario(gsOption, phsOption); + %% 1. define Constants Constants = io.define_constants(); - ParaPlant = io.define_plant_constants(SiteProperties, phwsfMethod); + RTB = 1000; % initial root total biomass (g m-2) if strcmp(SiteProperties.sitename, 'CH-HTC') @@ -103,10 +99,6 @@ % [Rl, Ztot] = Initial_root_biomass(RTB, ModelSettings.DeltZ_R, ModelSettings.rroot, ModelSettings.ML, SiteProperties.landcoverClass(1)); numSoilLayer = ModelSettings.ML; soilThickness = ModelSettings.DeltZ_R'; - [RootProperties,soilDepth] = calRootProperties(SiteProperties, ParaPlant, numSoilLayer, soilThickness, RTB); - Rl = RootProperties.lengthDensity; - Ztot = soilDepth; - soilDepthB2T = flipud(soilDepth); %% 2. simulation options @@ -116,6 +108,20 @@ parameter_file = {'input_data.xlsx'}; options = io.readStructFromExcel([path_input char(parameter_file)], 'options', 2, 1); + % Set senarios + gsOption = options.gsOption; + phsOption = options.plantHydraulics; + [gsMethod, phwsfMethod] = setScenario(gsOption, phsOption); + + % define parameters for plant hydraulics module + ParaPlant = io.define_plant_constants(SiteProperties, phwsfMethod); + + % calculate root parameters + [RootProperties,soilDepth] = calRootProperties(SiteProperties, ParaPlant, numSoilLayer, soilThickness, RTB); + Rl = RootProperties.lengthDensity; + Ztot = soilDepth; + soilDepthB2T = flipud(soilDepth); + if options.simulation > 2 || options.simulation < 0 fprintf('\n simulation option should be between 0 and 2 \r'); return diff --git a/src/setScenario.m b/src/setScenario.m index 6b831b9f..a08cfca1 100644 --- a/src/setScenario.m +++ b/src/setScenario.m @@ -4,11 +4,12 @@ hydraulic pathway option. input: - gsOption: a string indicate which stomatal conductance will be used + gsOption: a number indicate the stomatal conductance approach: + 1 for BallBerry, ref[1] + 2 for Medlyn, ref[2] phsOption: a number indicate whether use plant hydraulic pathway output: - biochemical: a string, indicate which function will be used gsMethod: a number indicate the stomatal conductance scheme: 1 for BallBerry 2 for Medlyn @@ -17,15 +18,25 @@ 2 for ED2, 3 for PHS + references: + [1] Ball, J. T., et al. (1987). A Model Predicting Stomatal Conductance + and its Contribution to the Control of Photosynthesis under + Different Environmental Conditions, Springer, Dordrecht. + + [2] Medlyn, B. E., et al. (2011). "Reconciling the optimal and empirical + approaches to modelling stomatal conductance." + Global Change Biology 17(6): 2134-2144. + %} %% Set scenario - %biochemical = @biochemical; - if strcmp(gsOption, 'BallBerry') - gsMethod = 1; - elseif strcmp(gsOption, 'Medlyn') - gsMethod = 2; + switch gsOption + case 1 + gsMethod = 1; % BallBerry stomatal conductance + case 2 + gsMethod = 2; % Medlyn stomatal conductance end + %% PHS setting if phsOption % if phsOption is true, PHS open. Set plant water stress method. switch phsOption