noaa-owp / evapotranspiration Goto Github PK

One of the standalone scripts, ./make_and_run_read_forcings.sh, includes a
Unit Test print out. Values do not match but expected values seem to be hard-coded? See pet.c here.

Reminder to look into this as potential bug or remove printout.

aerodynamic method:
calculated instantaneous potential evapotranspiration (PET) =8.977490e-08 m/s

 #----------- BEGIN UNIT TESTS   ---------------# 

 #-----------       UNIT TEST    ---------------# 
solar elevation angle is 47.376995 degrees,
 and should be: 43.29101185 degrees 

 #-----------       UNIT TEST    ---------------# 
solar azimuth angle is 213.491150 degrees,
 and should be: 225.06371958 degrees 

 #-----------       UNIT TEST    ---------------# 
solar local hour angle is 21.942448 degrees,
 and should be: 31.01549773 degrees 

 #-----------       UNIT TEST    ---------------# 
potential ET is 8.977490e-08 m/s,
 and should be: 8.977490e-08 m/s 

 #------------  END UNIT TESTS   ---------------# 

Is it worth considering renaming this repo to 'potential_ET' (or similar) to avoid confusion with (actual) evapotranspiration flux calculations, which occur in the hydrology model modules?

Based off get_time_step() def here, as well as any general search for time_step_size_s, seems like we are hard-coding an assumption for time step to be of units s.

Is this is case? If so, is this BMI compliant?

Reminder to further look into spec/implementation

A member of the pet_param struct, heat_transfer_roughness_length_m is NOT initialized in the BMI code. The related momentum_transfer_roughness_length_m is initialized in the BMI code, but BOTH of these parameters are used to determine (this condtional branch)[https://github.com/NOAA-OWP/evapotranspiration/blob/1df115fdd5376fe37be8ceafcb5df3b2dd0577e4/include/pet_tools.h#L408]. This value needs some initial value (0 based on the code in the pet_tools.h) before used in this function!

PET option 5, the Penman-Monteith method is optimized for relatively even grass/crop surfaces. In these areas, the zero plane displacement height, d, is typically less than the wind speed and/or relative humidity measurement height. This is important because estimating surface resistance, rs, requires taking the natural log of the measurement height, z, minus d. The natural log of a negative number (i.e., d > z) is NaN.

Equation from FAO.

Current behavior

When PET is run with option 5 in forested locations (d > z), it returns NaN. Part of this is because we are using AORC forcing data, which assumes z for wind is 10 m and z for humidity is 2 m. If running option 5 with above-canopy met towers, this would not be a problem.

Expected behavior

We should be able to estimate PET in forested locations using P-M as other researchers have done. We need to either modify the z estimate to be less than the measurement height or substitute empirically derived resistance values. E.g., 43–67 s/m from Amatya et al. (2016).

From Fred's comments in the Readme:
the skin temperature of the earth's surface, TODO: should come from another module to calc. soil or veg. temp.

More specifically, get_end_time() should return FLT_MAX when num_timesteps is not directly defined (or is 1). update_until() fails because def considers end_time as stopping condition.

See https://bmi.readthedocs.io/en/latest/#get-end-time

From Fred's comment in the Readme:

the incoming solar (shortwave) radiation. If not provided it is computed from d,e,f, using an updated method similar to the one presented in Bras, R.L. Hydrology. Requires value of the Linke atmospheric turbidity factor, which varies from 2 for clear mountain air to 5 for smoggy air. According to Hove & Manyumbu 2012, who calculated values over Zimbabwe that varied from 2.14 to 3.71. Other values exist in the literature. TODO: This turbidity factor could be calculated from satellite obs. or maybe NOAA already does this? All radiation calculations needed for 1, 3, 4, and 5 require net radiation calculations at the land surface. The net radiation is calculated using a, c, d, e, f, g, j, plus the Linke turbidity factor, which can be estimate from satellite observations.
For reference, from Readme:
a) temperature and (relative-humidity or specific humidity) and the heights at which they are measured.
b) near surface wind speed measurement and the height at which it was measured.
c) the ambient atmospheric temperature lapse rate
d) the fraction of the sky covered by clouds
e) (optional) the height above ground to the cloud base. If not provided, then assumed.
f) the day of the year (1-366) and time of day (UTC only!)
g) the skin temperature of the earth's surface.
h) the zero-plane roughness height of the atmospheric boundary layer assuming log-law behavior (from land cover)
i) the average root zone soil temperature, or near-surface water temperature in the case of lake evaporation.
j) the incoming solar (shortwave) radiation. If not provided it is computed from d,e,f, using an updated method similar to the one presented in Bras, R.L. Hydrology. Requires value of the Linke atmospheric turbidity factor, which varies from 2 for clear mountain air to 5 for smoggy air. According to Hove & Manyumbu 2012, who calculated values over Zimbabwe that varied from 2.14 to 3.71. Other values exist in the literature.
NOTE THE VALUE OF evapotranspiration_params.zero_plane_displacement_height COMES FROM LAND COVER DATA.

Current behavior

We currently house AORC forcing code in both the PET and CFE repos. This complicates maintenance and development as it needs to be performed twice.

Expected behavior

PET and CFE should each have an external submodule that points to the newly created aorc_bmi repo.

Steps to remedy behavior

1. First complete and close issue 1 on aorc_bmi to have the forcing code and BMI wrapper housed within the new repo
2. Fork the PET repo and start a new development branch
3. Remove the /forcing_code directory from the main level of the branch
4. Create a new /extern/forcing_code directory and subdirectory from the main level
5. Add the aorc_bmi submodule: git submodule add https://github.com/NOAA-OWP/aorc_bmi.git extern/forcing_code
6. Add a README to /extern
7. Modify the build and run scripts as necessary in PET
8. Test
9. Commit changes, push, and submit PR

Current behavior

PET calls for precipitation as part of its BMI specification, which is not necessary because the model does not use precipitation.

Expected behavior

Remove precipitation from input_var functions.

The time_step_size_s parameter is used to set the BMI model time step, which is used in Update_until to determine how many Update steps to call. If this parameter isn't provided in the configuration, it is left uninitialized and leads to undefined behavior. For example, if the compiler assigns 0 for the uninitialized variable, then in Update_until there is a divide by 0, which can cause a crash, or worse, can make n_steps become inf (again, this is undefined behavior...) which causes an infinite loop.

This either needs to be a required parameter and validated to not be 0, or an optional value (also validated to not be 0) and a non-zero default assigned during initialization.

PET options (2) aerodynamic, (3) combined, and (5) Penman-Monteith methods are currently providing unreasonable PET values.

Current behavior

Table below shows PET mean annual values (mm/year) for selected CAMELS basins for the five PET options. PET options 2, 3 and 5 present unrealistically high PET values.

Possible causes: unit issues for the input or for the parameters in the config file.

Short description explaining the high-level reason for the new issue.
Some forcing variables get passed in through BMI incorrectly

Current behavior

examples:
specific_humidity_2m_kg_per_kg 36893488147419103232.000000
LWDOWN after set value -36893488147419103232.000000
u_wind_speed_10m_m_per_s -36893488147419103232.000000
v_wind_speed_10m_m_per_s -36893488147419103232.000000

Expected behavior

realistic forcing values get passed in through BMI

Steps to replicate behavior (include URLs)

1. check the variable types, make sure doubles pass to doubles.

Looking at the provided example configuration files and the readme instructions for calculating zero plane roughness length, the only way i'm coming up with something close to 0.0003 is by linearly interpolating between the values in the roughness class roughness length(m) lookup table. Is that how the look up table is supposed to be used?

vegetation_height_m=0.12
zero_plane_displacement_height_m=0.0003

source

zero_plane_displacement_height_m = y1 + (x - x1) * (y2 - y1) / (x2 - x1)
                                 = 0.0002 + (0.12 - 0) * (0.0005 - 0.0002) / (0.2 - 0)
                                 = 0.00038

From an admittedly pedantic reading of the readme instructions for calculating zero plane roughness, the word selected did not tip me off interpolate between the table class values.

The zero plane roughness length,"d" can be approximated as 2/3 of the vegetation height (H): d=2/3*H for grassland/cropland or selected from above table based on comprehensive land cover categories.

source