feat: adding LandSurfaceTemperature module for LST modeling #633
Conversation
The module sets up the first keyword parameter. serves as foundation to implement the methods from the Devi Purnamasari et al. 2025 doi: https://doi.org/10.5194/hess-29-1483-2025
-add read_lst_inputs function handling global data ie radiation, albedo and crop type -probably add emissivity later
- read inputs from read lst inputs - initialize and empty grid
…read in init func
-requires {model:{lst__flag:true}}
-following SnowModel if do_lst will initialize LSTModel or not
Add energy balance calculations with radiation and heat fluxes Integrate atmospheric forcing and SBM soil model coupling Implement detailed LST calculations following Purnamasari et al. 2025
vers-w
left a comment
vers-w
left a comment
There was a problem hiding this comment.
Nice start with this new feature @michaelohanrahan !
I added some comments, main suggestion is to include this new feature as part of LandHydrologySBM, make use of variables and parameter fields in the model struct, to move some parameters or forcing to other structs (e.g. for shared land parameters) and to check/discuss what hydromt_wflow should compute.
Add albedo, shortwave radiation (RS_in), and 2m wind speed (u2m) to forcing Support land surface temperature model requirements
Add LST, radiation fluxes, and aerodynamic resistance with proper units Maintain Celsius temperature scale for land surface temperature
…d surface temperature - Use explicit LandSurfaceTemperature* names for model, parameters, and variables - Clean up struct definitions for clarity and maintainability - Update to handle expanded atmospheric forcing structure
Add LST model update call in main timestep loop after soil model updates
Include land_surface_temperature flag in model configuration for SBM GWF model
Include land_surface_temperature flag in model configuration for SBM model
vers-w
left a comment
vers-w
left a comment
There was a problem hiding this comment.
Thanks for the changes @michaelohanrahan ! Some additional comments from my side.
- Remove LandSurfaceTemperatureParameters struct entirely - Update LandSurfaceTemperatureModel to only contain variables (no parameters) - Modify update_timestep_land_surface_temperature function to accept network::NetworkLand and vegetation_parameters::VegetationParameters as parameters - Update function to use network.latitude[i] instead of land_surface_temperature_model.parameters.latitude[i] - Update function to use vegetation_parameters.canopy_height[i] instead of land_surface_temperature_model.parameters.crop_height[i] - Add aerodynamic resistance calculation within the update function - Fix bug in calculate_land_surface_temperature function to use correct constant name cp - Update both update! function signatures to include new parameters
- Update all land surface temperature variable mappings to use new variable names: * LST → land_surface_temperature * RSN → net_shortwave_radiation * RLN → net_longwave_radiation * Rnet → net_radiation * LE → latent_heat_flux * H → sensible_heat_flux * ra → aerodynamic_resistance
…on params. Removing clock and dt from model update
…opy types - Implemented a robust aerodynamic resistance and wind speed correction function following Allen et al. (1998) and FAO56. - The new approach uses the standard log-profile for short canopies (e.g., grass, wheat) when the wind measurement height is well above the displacement height. - For tall canopies (e.g., forests, tall crops), the function now caps the displacement height (d ≤ 0.9 * z_measured) and uses lower d/h ratios (0.5–0.6) to avoid invalid log-profile math. - If the log-profile is not valid (e.g., measurement height below or near d), the function falls back to using the measured wind speed and an empirical resistance (e.g., ra = 208/u2). - Added warnings and debug output to inform users when fallbacks are used. - This change ensures physical consistency and prevents negative/undefined aerodynamic resistance values for all crop and canopy types.
… empty nan array.
…nopy and also units are correctly in range.
Co-authored-by: Willem van Verseveld <[email protected]>
…into landsurfacetemp
| # Read canopy height | ||
| lens = lens_input_parameter(config, "vegetation__canopy_height") | ||
| canopy_height = | ||
| ncread(dataset, config, lens; sel = indices, defaults = 0.12, type = Float64) |
There was a problem hiding this comment.
Any reason to use a value of 0.12 here? The default value as part of the Sediment model is 0.5.
There was a problem hiding this comment.
yes i use this to default to the FAO (Allen '98) grass reference as default
| dh_ratio = 2.0 / 3.0 | ||
| z0h_ratio = 0.2 | ||
| else | ||
| error("Canopy height $canopy_height is not supported") |
There was a problem hiding this comment.
Wondering if this error message is really needed?
What height would not be supported as it is already set to a minimum of 0.12? If this function is only called for cells with the correct range of canopy_height (see also comment above) I think catching this error is not needed. Raising an error in a function that is called in a for loop and for each timestep can get computationally expensive and also for a user it may get overwhelming, depending how many times the error is raised.
There was a problem hiding this comment.
Agreed! A remnant of testing that should now be removed. It was indeed an overwhelming error at one point, so I completely understand the suggestion.
|
|
||
| # Ensure positive aerodynamic resistance | ||
| if ra <= 0 | ||
| error("Aerodynamic resistance is negative: $ra") |
There was a problem hiding this comment.
Instead of throwing an error I would recommend to solve this with a minimum value for ra or if possible resolve this before ra is computed.
There was a problem hiding this comment.
I don't think it is possible any longer to calculate a negative with the floor in wind speed, and the handling of the taller canopies. I would be hesitant to assuming a defualt value for ra as the end result is quite sensitive and the combination of vars that lead to negative values are not so predictable. I will retain for testing and add a TODO to remove after I am confident the issue is eliminated.
| if isnan(aerodynamic_resistance) || | ||
| isinf(aerodynamic_resistance) || | ||
| aerodynamic_resistance <= 0 | ||
| @warn "Invalid aerodynamic resistance: $aerodynamic_resistance, using air temperature as land surface temperature" |
There was a problem hiding this comment.
As mentioned above probably better to leave out log messages in functions that are applied in a for loop for each cell and at each timestep. Would recommend to give it a minimum value (or not including cells that are invalid for these kind of computations).
There was a problem hiding this comment.
is it allowed to return a nan (or nodata) value in the case of an invalid iteration?
…into landsurfacetemp
committing suggestion Co-authored-by: Willem van Verseveld <[email protected]>
Explanation
Working towards implementation of the methods from the Devi Purnamasari et al. 2025 paper: "Identifying irrigated areas using land surface temperature and hydrological modelling: application to the Rhine basin" doi
Implementation details:
Checklist
masterAdditional Notes (optional)
Add any additional notes or information that may be helpful.