<p>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!</p>
<p>Original question:</p>
<p>I am trying to model Penman-Monteith reference evapotranspiration in MATLAB using Ameriflux data (<a href="https://ameriflux.lbl.gov/sites/siteinfo/US-Ro4" rel="nofollow noreferrer">https://ameriflux.lbl.gov/sites/siteinfo/US-Ro4</a>, 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.</p>
<p>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.</p>
<p>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.</p>
<pre><code>clear
load Project2_AmerifluxData_2021.mat
dates = datevec(FluxData.TIMESTAMP_START);
selector = dates(:, 1) == 2021 &amp; dates(:, 2) == 6 &amp; 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 &amp; 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
</code></pre>
<p><a href="https://i.sstatic.net/9a4bC.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/9a4bC.png" alt="Table of data variable descriptions" /></a>
<a href="https://i.sstatic.net/jaeNm.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/jaeNm.png" alt="Data variables in MATLAB" /></a>
<a href="https://i.sstatic.net/1AQzr.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/1AQzr.png" alt="More data variables in MATLAB" /></a>
<a href="https://i.sstatic.net/XXNrt.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/XXNrt.png" alt="Graph example 1" /></a>
<a href="https://i.sstatic.net/a7yga.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/a7yga.png" alt="Graph example 2" /></a>
<a href="https://i.sstatic.net/fnp6k.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/fnp6k.png" alt="Graph example 3" /></a></p>

CC BY-SA 4.0
2024-03-13T08:15:35.693Z
Hi,pls ask on other stack exchange Pages (e.g. stack overflow) I Think you will finde there faster an answer :)
Weiss
CC BY-SA 4.0
2024-03-13T08:42:19.050Z
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!
Stan Nonna
CC BY-SA 4.0
2024-03-13T08:58:08.520Z
If you ask at the outher Page pls deleted the Questionen Here. Outherwise they wount Help you
Weiss
CC BY-SA 4.0
2024-03-13T09:18:10.003Z
OK, I might keep the question here instead because I don't see a hydrology tag in stack overflow.
Stan Nonna
CC BY-SA 4.0
2024-03-13T09:20:44.680Z
But to bei honest i dont think you will get an answer Here soon. I personaly would try IT at Stack overflow and If i wouldnt get an answer in lets say 3 days i would cone Back :)
Weiss
CC BY-SA 4.0
2024-03-13T09:21:38.897Z
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.
Jean-Marie Prival
CC BY-SA 4.0
2024-03-13T09:29:25.143Z
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 :)
JeopardyTempest
CC BY-SA 4.0
2024-03-13T09:40:00.417Z
I won't be able to debug the code and make a proper answer... But just having a quick look at it, some things seem odd. Like when you compute your Delta, you comment that its unit is Pa/K, but you use Flux_6_17_21.TA, which is in °C, instead of TA_K. Is it normal?
Jean-Marie Prival
CC BY-SA 4.0
2024-03-13T09:47:56.003Z
I also don't think SO will get very far for this question.  It'll probably drown in the huge number of posts there and go unnoticed.
gerrit
CC BY-SA 4.0
2024-03-13T09:58:12.890Z
@Jean-MariePrival Thank you! I have done my best to check the units--I agree that they're usually the culprit when something goes wrong. I still could have missed something. I believe delta is OK. That parameter is the rate of change of saturation-specific humidity with air temperature. Because the difference between the starting and ending temperature is the same regardless of whether temp is in K or °C, I didn't think it would be an issue to switch units at that point.
Stan Nonna
CC BY-SA 4.0
2024-03-24T12:30:28.287Z
If you were able to solve this yourself in the meantime, maybe post and accept your answer so that this does not appear unanswered? I was just about to deep-dive into PM when I noticed your "update" in the question.
dimfalk
<hydrology><models><evaporation><matlab>
What's wrong with my Penman-Monteith model?
Stan Nonna
# What's wrong with my Penman-Monteith model?

(Asked by: username_0 on 2024-03-13)

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
load Project2_AmerifluxData_2021.mat
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

```

[![Table of data variable descriptions](https://i.sstatic.net/9a4bC.png)](https://i.sstatic.net/9a4bC.png)
[![Data variables in MATLAB](https://i.sstatic.net/jaeNm.png)](https://i.sstatic.net/jaeNm.png)
[![More data variables in MATLAB](https://i.sstatic.net/1AQzr.png)](https://i.sstatic.net/1AQzr.png)
[![Graph example 1](https://i.sstatic.net/XXNrt.png)](https://i.sstatic.net/XXNrt.png)
[![Graph example 2](https://i.sstatic.net/a7yga.png)](https://i.sstatic.net/a7yga.png)
[![Graph example 3](https://i.sstatic.net/fnp6k.png)](https://i.sstatic.net/fnp6k.png)

**Tags:** hydrology, models, evaporation, matlab

> **Comments:**
> **username_1**: Hi,pls ask on other stack exchange Pages (e.g. stack overflow) I Think you will finde there faster an answer :)  
> **username_0**: 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!  
> **username_1**: If you ask at the outher Page pls deleted the Questionen Here. Outherwise they wount Help you  
> **username_0**: OK, I might keep the question here instead because I don't see a hydrology tag in stack overflow.  
> **username_1**: But to bei honest i dont think you will get an answer Here soon. I personaly would try IT at Stack overflow and If i wouldnt get an answer in lets say 3 days i would cone Back :)  
> **username_2**: 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.  
> **username_3**: 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 username_2 is probably your best hope :)  
> **username_2**: I won't be able to debug the code and make a proper answer... But just having a quick look at it, some things seem odd. Like when you compute your Delta, you comment that its unit is Pa/K, but you use Flux_6_17_21.TA, which is in °C, instead of TA_K. Is it normal?  
> **username_4**: I also don't think SO will get very far for this question.  It'll probably drown in the huge number of posts there and go unnoticed.  
> **username_0**: @username_2 Thank you! I have done my best to check the units--I agree that they're usually the culprit when something goes wrong. I still could have missed something. I believe delta is OK. That parameter is the rate of change of saturation-specific humidity with air temperature. Because the difference between the starting and ending temperature is the same regardless of whether temp is in K or °C, I didn't think it would be an issue to switch units at that point.  
> **username_5**: If you were able to solve this yourself in the meantime, maybe post and accept your answer so that this does not appear unanswered? I was just about to deep-dive into PM when I noticed your "update" in the question.  

---

## 0 answers

