# What's wrong with my Penman-Monteith model?

UPDATE: I believe I have found the mistake. I had missed that the project instructions stated that conductance of sensible heat flux should be equal to aerodynamic conductance. The model behaves as expected now that I've made that adjustment. Thank you to everyone who chimed in!

Original question:

I am trying to model Penman-Monteith reference evapotranspiration in MATLAB using Ameriflux data (https://ameriflux.lbl.gov/sites/siteinfo/US-Ro4, data from June 17, 2021). The task is to calibrate the model to the evapotranspiration data (the latter calculated using latent heat flux) by adjusting the value for canopy conductance. The canopy conductance was expected to lie somewhere between 0.002 and 0.03 m/s. I was asked to plot an overestimate, an underestimate, and a best fit model on the same plot as the ET data.

Unfortunately, it seems that the model is off. Regardless of the choice of G_c, the plot oscillates in unexpected places, and contains negative values. Also, the G_c value that minimizes error--as far as I've found with the current approach of manual adjustment--lies outside the recommended range.

I will owe a life debt to anyone who can point me toward my mistake(s). Apologies for the needlessly bothersome aspects of the code, the large images, and any breaches I may have made of forum etiquette: I'm a newbie to both MATLAB and StackExchange.

clear
dates = datevec(FluxData.TIMESTAMP_START);
selector = dates(:, 1) == 2021 & dates(:, 2) == 6 & dates(:, 3) == 17;
Flux_6_17_21 = FluxData(selector, :);

% Calculate evapotranspiration from LE data
Lv = 2257000; % Latent heat of vaporization [J/kg]
E = Flux_6_17_21.LE/Lv; % Evapotranspiration rate from LE data [kg/(s·m^2)]

% Account for other terms in the Penman-Monteith equation

% Temperature
TA_K = Flux_6_17_21.TA + 273.15; % Air temperature [K] at measurement height z
TS_all = [Flux_6_17_21.TS_1_1_1, Flux_6_17_21.TS_1_2_1, Flux_6_17_21.TS_1_3_1, Flux_6_17_21.TS_1_4_1];
TS_C = mean(TS_all, 2); % Mean of different sensors' surface temperature readings [deg C]
TS_K = TS_C + 273.15; % Mean surface temperature [K]

% Vapor pressure
PA = 1000 * Flux_6_17_21.PA; % Atmospheric pressure [Pa]
a = 611; % Clausius-Clapeyron fitting parameter [Pa]
b = 17.27; % C-C parameter [-]
c = 237.3; % C-C parameter [deg C]
e_ss = a * exp(b * TS_C ./ (TS_C + c)); % Saturated water vapor
% pressure [Pa] at surface, from Clausius-Clapeyron
e_sz = a * exp(b * Flux_6_17_21.TA ./ (Flux_6_17_21.TA + c)); % Sat
% water vapor pressure [Pa] at z
e_az = e_sz .* Flux_6_17_21.RH ./ 100; % Actual water vapor pressure [Pa] at z
D = Flux_6_17_21.VPD_PI .* 100; % Vapor pressure deficit [Pa]

% Delta [Pa/K]
delt = a * b * c ./ (Flux_6_17_21.TA + c).^2 .* exp(b * Flux_6_17_21.TA ./ (Flux_6_17_21.TA + c));
% Slope of C-C eqn

% Air density
R = 8.314; % Universal gas constant [J/(K·mol)]
c_pd = 1005; % Specific heat capacity of dry air [J/(kg·K)]
c_pw = 1930; % Specific heat capacity of water vapor [J/(kg·K)]
mol_fract = e_az ./ PA; % Mol fraction [-] of water vapor at z
c_p = (1-mol_fract) .* c_pd + mol_fract .* c_pw; % Specific heat capacity
% of moist air [J/(kg·K)]
p_d = PA - e_az; % Partial pressure of dry air [Pa] at z
M_a = 0.02897; % Molar mass of dry air [kg/mol]
M_w = 0.018; % Molar mass of water [kg/mol]
rho_a = (M_w * e_az + M_a * p_d) ./ (R * TA_K); % Mass density of air [kg/m^3]

% Aerodynamic conductance
kappa = 0.41; % von Karman constant [-]
z = 10; % Measurement height [m]
h = 0.5; % Height of vegetation canopy [m]
d = 0.667 * h; % Zero plane displacement [m]
zm = 0.123 * h; % Roughness length for momentum transfer [m]
zh = 0.1 * zm; % Roughness length for heat & water vapor [m]
G_a = kappa^2 * Flux_6_17_21.WS ./ (log((z - d)/zm) * log((z - d)/zh)); % Aerodynamic
% conductance [m/s]

% Other conductances (all [m/s])
G_c = 0.12; % Canopy conductance, used interchangeably with surface conductance. Vary this to find best fit
G_L = G_c .* G_a ./ (G_c + G_a); % Total conductance of latent heat flux
G_H = Flux_6_17_21.H ./ (rho_a .* c_p .* (TS_K - TA_K)); % Conductance of sensible heat flux

% Psychrometric constant [Pa·K]
gamma = c_p .* PA / (Lv * M_w / M_a);

% Recalculate ET with Penman-Monteith equation
E_PM = (delt .* (Flux_6_17_21.NETRAD - Flux_6_17_21.G) + rho_a .* c_p .* G_H .* D) ./ (Lv ...
* (delt + gamma .* G_H ./ G_L)); % [kg/(s·m^2)]

Error = mean(rmse(E, E_PM, 2)) % Look for value of G_c that minimizes error between E and E_PM

% Variation #1
G_c_1 = 0.03; % [m/s]
G_L_ov = G_c_1 .* G_a ./ (G_c_1 + G_a); % [m/s]
E_PM_ov = (delt .* (Flux_6_17_21.NETRAD - Flux_6_17_21.G) + rho_a .* c_p .* G_H .* D) ./ (Lv ...
* (delt + gamma .* G_H ./ G_L_ov)); % [kg/(s·m^2)]

% Variation #2
G_c_2 = 0.002; % [m/s]
G_L_und = G_c_2 .* G_a ./ (G_c_2 + G_a); % [m/s]
E_PM_und = (delt .* (Flux_6_17_21.NETRAD - Flux_6_17_21.G) + rho_a .* c_p .* G_H .* D) ./ (Lv ...
* (delt + gamma .* G_H ./ G_L_und)); % [kg/(s·m^2)]

figure(1)
plot (Flux_6_17_21.TIMESTAMP_START, E, '-b')
hold on
plot (Flux_6_17_21.TIMESTAMP_START, E_PM, '--g')
plot (Flux_6_17_21.TIMESTAMP_START, E_PM_ov, '-.k')
plot (Flux_6_17_21.TIMESTAMP_START, E_PM_und, ':r')
datetick('x', 'HHPM')
xlabel('Hour')
ylabel('Evapotranspiration rate [kg/(s·m^{2})]')
legend ('Data', ['Penman-Monteith equation with G_c = ', num2str(G_c), ...
' m/s'], ['Penman-Monteith equation with G_c = ', num2str(G_c_1), ...
' m/s'], ['Penman-Monteith equation with G_c = ', num2str(G_c_2), ' m/s'])
legend('Location', 'northoutside')
hold off


• Hi,pls ask on other stack exchange Pages (e.g. stack overflow) I Think you will finde there faster an answer :) Commented Mar 13 at 8:15
• Oh, thank you for the tip! I wasn't aware of stack overflow. I don't know if I've messed up due to a misunderstanding of MATLAB commands or of the Penman-Monteith equation, or both. If it's the former, hopefully the other site can help! Commented Mar 13 at 8:42
• OK, I might keep the question here instead because I don't see a hydrology tag in stack overflow. Commented Mar 13 at 9:18
• I'd say: begin by double and triple checking your units. You do a lot of conversions, for °C to K, from hPa or kPa to Pa... From my experience, when writing MATLAB script like this, 95% of the time errors will come from a unit conversion issue. Commented Mar 13 at 9:21
• Generally it's not suggested to bounce a question around and try a lot of places. I agree that the limited community here will probably(?) take a while to respond if at all. Then again, not sure SO will have the knowledge base to help. An answer from someone like Jean-Marie is probably your best hope :) Commented Mar 13 at 9:29