Spatially varying canopy height #1502

braghiere · 2025-10-22T17:59:25Z

Purpose

Add support for spatially-varying canopy height in ClimaLand vegetation model, replacing the fixed 1m default with realistic height data from CLM artifacts (capped at 8m).

Content

Core Implementation

Canopy height parameterization: Added clm_canopy_height() function to read spatially-varying canopy height from CLM artifact data (src/standalone/Vegetation/spatially_varying_parameters.jl)
Height capping: Implemented effective_canopy_height() function to cap heights at 8m for numerical stability while preserving spatial patterns
Type system updates: Extended PrescribedBiomassModel to support spatially-varying height through type parameter HTH (canopy height type)
Boundary layer integration: Updated displacement height and roughness length calculations to properly broadcast with spatially-varying height Fields (src/standalone/Vegetation/canopy_boundary_fluxes.jl)
Bug fixes:
- Fixed h_leaf synchronization with canopy_height in PlantHydraulics
- Resolved Field comparison performance issues in height assertions
- Fixed integer overflow in effective_canopy_height Boolean operations
- Corrected GPU Boolean Field handling

Testing and Validation

Unit tests: Comprehensive test suite covering spatially-varying height functionality, including Field operations, capping behavior, and integration with canopy model

Impact Assessment

Figure: Sensible heat flux comparison showing Run 1 (default 1m height), Run 2 (spatially-varying height), absolute difference, and percent change.

The spatially-varying canopy height primarily affects:

Turbulent fluxes: Changes in aerodynamic resistance modify sensible and latent heat exchange with atmosphere
Regional patterns: Strongest impacts in forested regions (tropical/temperate) where height exceeds 1m default
Minimal impact: Grasslands and croplands show small changes (CLM heights ~1m already)

To-do

Add NEWS.md entry after reviewer feedback

I have read and checked the items on the review checklist.

Previously h_leaf was hardcoded to FT(1), now it uses toml_dict["canopy_height"] to ensure consistency between biomass and hydraulics height parameters.

…acts - Reads spatially varying canopy height from CLM vegetation_properties_map.nc - Returns ClimaCore Field with height values for each grid cell - Follows same pattern as clm_rooting_depth() and clm_medlyn_g1() - Heights are based on dominant PFT from CLM5 MONTHLY_HEIGHT_TOP data

CLM data uses 'z_top' (Canopy top height) not 'canopy_height'. Data investigation findings: - 192x288 grid (0.9x1.25 degrees) - Height range: 0-40.13m - 92.5% < 1m (grasslands/crops) - 7.4% >= 10m (4068 cells, tall forests) - Max height: 40.13m (tropical forests) - Median: 0m (sparse/no vegetation dominant)

- Add HTH type parameter to support Union{FT, Field} for height - Update struct definition in biomass.jl to accept both scalar and Field height - Update outer constructor to pass height type to struct instantiation - Update Canopy.jl constructor to include height type parameter - This allows height to be either a scalar (backward compatible) or a Field (new functionality)

- Implement effective_canopy_height() to cap heights below atmospheric reference - Default buffer of 2m ensures adequate space for flux calculations - Logs warning with statistics when capping occurs - Addresses constraint that ERA5 forcing is at 10m reference height - Prevents numerical issues from canopy extending above atmospheric layer

- Import clm_canopy_height and effective_canopy_height functions - Read canopy height from CLM artifact data in snowy_land_pmodel.jl - Apply capping to 8m (10m atmospheric height - 2m buffer) - Pass spatially-varying height to PrescribedBiomassModel - This replaces the previous fixed 1m canopy height with realistic CLM values

Tests added: 1. spatial_parameters.jl: - Test clm_canopy_height() reads CLM data correctly - Test effective_canopy_height() capping with default and custom buffers - Verify proper field structure and value constraints 2. test_spatially_varying_canopy_height.jl (new file): - Test PrescribedBiomassModel with scalar height (backward compatibility) - Test PrescribedBiomassModel with Field height (new functionality) - Test integration with CLM height data and capping - Test full CanopyModel construction with spatially-varying height All tests verify: - Correct Field axes and types - Height constraints (capping below atmospheric reference) - Non-negative heights - Backward compatibility with scalar heights

Create experiment script to compare: - Run 1: Default configuration with height=1m everywhere - Run 2: Spatially-varying height from CLM data (capped at 8m) Both runs use same 1-year period (2008-03-01 to 2009-03-01) to enable direct comparison of model outputs and assess impact of spatially-varying canopy height on simulated fluxes and states.

- Replace count() with ifelse() to avoid Boolean Fields on GPU - Use float mask (1.0/0.0) instead of boolean for counting - Fixes 'Bool cannot be represented using Float64' error on CUDA

- Remove Int() conversion that caused overflow with large GPU arrays - Keep n_capped_sum as Float for display - Use round() for cleaner output

…t_height() - Changed FT(0.67) * height to FT(0.67) .* height - Fixes MethodError when height is a Field (spatially-varying) - Maintains backward compatibility with scalar height

- Changed biomass.height == hydraulics_height to use all(isapprox.(...)) - Original == with Field was creating Field of booleans (very slow) - Maintains backward compatibility with scalar height using ≈ - This was causing 3x+ slowdown in spatially-varying height runs

These simulation outputs have been moved to ~/climaland_outputs/ for local analysis and are not part of the core feature implementation.

Wrap parent() array accesses with Array() to convert from GPU to CPU arrays before scalar indexing. This prevents 'Scalar indexing is disallowed' errors when running tests on GPU.

kmdeck · 2025-10-22T20:26:43Z

src/standalone/Vegetation/canopy_boundary_fluxes.jl

 """
 function ClimaLand.displacement_height(model::CanopyModel{FT}, Y, p) where {FT}
-    return FT(0.67) * model.biomass.height
+    return FT(0.67) .* model.biomass.height


Any operations like:
a .* b will allocate, so we should be careful to not do this in functions that are evaluated each timestep (it is OK if we do it selectively in initialization).

Here maybe we can do:
@. lazy( FT(0.67) * model.biomass.height)

which will only execute the operation when it is part of a larger set of operations updating something in place (no allocation if in place)

kmdeck · 2025-10-22T20:28:05Z

src/standalone/Vegetation/Canopy.jl

+                isapprox.(biomass.height, hydraulics_height, rtol = 1e-10),
+            ) "Biomass height Field must match hydraulics height"
+        else
+            @assert biomass.height ≈ hydraulics_height "Biomass height must match hydraulics height"


I dont quite understand how this works, because the hydraulics_height is a scalar. How can this pass when biomass.height is spatialy varying?

In general we need to go through src and make sure these are consistent, because in some places we use the 'hydraulics" height as the height

It would be good to document what biomass height means and what hydraulics height means, and what they're each used for

You're right - the hydraulics compartments are always scalar (1m). When biomass.height varies spatially, only the aerodynamic calculations use those varying heights. The hydraulics structure stays uniform everywhere.

This is a current limitation - we can extend hydraulics to support spatial variation in a follow-up PR.

kmdeck · 2025-10-22T20:29:21Z

experiments/long_runs/snowy_land_pmodel.jl

+        effective_canopy_height(raw_canopy_height, FT(10.0); buffer = FT(2.0))
+
+    # Create biomass model with spatially-varying height
+    biomass = PrescribedBiomassModel{FT}(


I think it would be good to add lines 112-116 into the PrescribedBiomassModel constructor, so it is not something that has to be updated each time (if we decide to include this as the default)

kmdeck

let's try running this and see how the SHF, LHF, etc change compared to observations

Moving forwards, I think there are two ways this can be extended:

make a canopy roughness model type. right now we are using physical quantities (height) to compute roughness and displacement height (using physical equations someone has come up with). Something more flexible would be to predict roughness, displacement height, from inputs like LAI, canopy height, etc using a data driven model. Having a type for "canopy roughness" would let us switch between the two easily. (I am separately adding this already in another PR)
in Biomass - create a type of "height model": this would be prescribed (constant or spatially varying) for now, PrescribedCanopyHeight. Then in the future we can also easily make different parameterizations for canopy height (e.g. prognostic, predicting it from input/time varying data). but this seems very far in the future compared to the first point

juliasloan25 · 2025-10-22T21:25:44Z

test/standalone/Vegetation/test_spatially_varying_canopy_height.jl

+    # Create a spatially-varying height field
+    coords = ClimaCore.Fields.coordinate_field(surface_space)
+    field_height = coords.lat .* 0 .+ FT(3.0)  # Constant field of 3m for testing


This is a little simpler, but what you have is fine too

Suggested change

# Create a spatially-varying height field

coords = ClimaCore.Fields.coordinate_field(surface_space)

field_height = coords.lat .* 0 .+ FT(3.0) # Constant field of 3m for testing

# Create a spatially-varying height field

field_height = zeros(surface_space) .+ FT(3) # Constant field of 3m for testing

juliasloan25 · 2025-10-22T21:31:28Z

src/standalone/Vegetation/Canopy.jl

+                isapprox.(biomass.height, hydraulics_height, rtol = 1e-10),
+            ) "Biomass height Field must match hydraulics height"
+        else
+            @assert biomass.height ≈ hydraulics_height "Biomass height must match hydraulics height"


It would be good to document what biomass height means and what hydraulics height means, and what they're each used for

juliasloan25 · 2025-10-22T21:49:32Z

src/standalone/Vegetation/spatially_varying_parameters.jl

+    canopy_height = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "vegetation_properties_map.nc"),
+        "z_top",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc, interpolation_method),
+    )


We could provide the effective_canopy_height as a preprocessing function instead, so we don't have to think about it wherever we access the canopy height data. This would be safer because we'll avoid users forgetting to apply the clipping.

file_reader_kwargs = (; preprocess_func = (data) -> effective_canopy_height(data, z_atm; buffer),),

And then we take z_atm and buffer as inputs to this function

effective_canopy_height would also have to change to take in a float canopy_height instead of the whole field, since here it'll be broadcasted over the field. We wouldn't be able to log the information you have in the warning, but I think that's okay. Usually when we add warnings like this we regret it later because it clutters the simulation output

juliasloan25 · 2025-10-22T21:54:47Z

experiments/long_runs/run2_spatially_varying_only.jl

@@ -0,0 +1,187 @@
+# RUN 2 ONLY: Spatially-Varying Canopy Height


If this is going to be merged into main it would be good to rename to something like snowy_land_varying_height.jl

juliasloan25 · 2025-10-22T21:56:57Z

experiments/long_runs/compare_canopy_heights.jl

+# module load climacommon
+# export CLIMACOMMS_DEVICE=CUDA
+# export CLIMACOMMS_CONTEXT=SINGLETON
+# cd /home/renatob/ClimaLand.jl


can we remove this line since it's specific to your directory?

juliasloan25 · 2025-10-22T22:00:05Z

test/standalone/Vegetation/spatial_parameters.jl

+    # Create a test field with known values including some that exceed the cap
+    coords = ClimaCore.Fields.coordinate_field(surface_space)
+    test_heights = coords.lat .* 0 .+ FT(5.0)  # Initialize to 5m


Same comment I left below -
This is a little simpler, but what you have is fine too

Suggested change

# Create a test field with known values including some that exceed the cap

coords = ClimaCore.Fields.coordinate_field(surface_space)

test_heights = coords.lat .* 0 .+ FT(5.0) # Initialize to 5m

# Create a spatially-varying height field

field_height = zeros(surface_space) .+ FT(5) # Initialize to 5m

juliasloan25 · 2025-10-22T22:02:59Z

test/standalone/Vegetation/test_spatially_varying_canopy_height.jl

@@ -0,0 +1,248 @@
+using Test


Thank you for writing such good tests!! This file should be added to test/runtests.jl so it runs with the test suite

…ng, clarify docs, update tests, and clean up experiment scripts

…ation

AlexisRenchon · 2025-10-29T19:50:02Z

Project.toml

 SurfaceFluxes = "0.12.3, 0.13"
 Thermodynamics = "0.14.1, 0.15"
 julia = "1.10"
+JuliaFormatter = "2.1.7"


I don't recall why exactly, but we still use JuliaFormatter 1. (not 2.)

JuliaFormatter shouldn't be a dependency of ClimaLand, you can keep it in your default environment

do:
]add JuliaFormatter@1

this will install the latest v1

to reset this to main, do:
git checkout origin/main Project.toml

braghiere · 2025-10-29T20:24:50Z

let's try running this and see how the SHF, LHF, etc change compared to observations

Moving forwards, I think there are two ways this can be extended:

make a canopy roughness model type. right now we are using physical quantities (height) to compute roughness and displacement height (using physical equations someone has come up with). Something more flexible would be to predict roughness, displacement height, from inputs like LAI, canopy height, etc using a data driven model. Having a type for "canopy roughness" would let us switch between the two easily. (I am separately adding this already in another PR)

in Biomass - create a type of "height model": this would be prescribed (constant or spatially varying) for now, PrescribedCanopyHeight. Then in the future we can also easily make different parameterizations for canopy height (e.g. prognostic, predicting it from input/time varying data). but this seems very far in the future compared to the first point

Largest changes are observed:

Sensible Heat (ERA5): Bias ↓ (8.27 → 7.89 W m⁻²) ~5% improvement, RMSE same.
here

Evaporative Fraction (ERA5): Bias slightly closer to 0 (−0.0626 → −0.0601) ~4% improvement, Phase same. (No RMSE.)
here

Surface Upward LW (CERES): Bias closer to 0 (−5.85 → −5.52 W m⁻²) ~6% improvement, RMSE ~ +0.1, Phase ≈ same.
here

All the other variables have minimal changes.

…ight

braghiere added 14 commits October 21, 2025 15:31

Fix: Sync PlantHydraulics h_leaf with biomass canopy_height

fd9893f

Previously h_leaf was hardcoded to FT(1), now it uses toml_dict["canopy_height"] to ensure consistency between biomass and hydraulics height parameters.

Fix GPU Boolean Field issue in effective_canopy_height

e3585e5

- Replace count() with ifelse() to avoid Boolean Fields on GPU - Use float mask (1.0/0.0) instead of boolean for counting - Fixes 'Bool cannot be represented using Float64' error on CUDA

Fix Int overflow in effective_canopy_height

147833a

- Remove Int() conversion that caused overflow with large GPU arrays - Keep n_capped_sum as Float for display - Use round() for cleaner output

Fix: broadcast scalar multiplication with height Field in displacemen…

e604889

…t_height() - Changed FT(0.67) * height to FT(0.67) .* height - Fixes MethodError when height is a Field (spatially-varying) - Maintains backward compatibility with scalar height

Clean up: Remove temporary log files and duplicate experiment scripts

1948ffc

Remove comparison experiment outputs from repository

79c760b

These simulation outputs have been moved to ~/climaland_outputs/ for local analysis and are not part of the core feature implementation.

braghiere added the enhancement New feature or request label Oct 22, 2025

braghiere requested review from AlexisRenchon and kmdeck October 22, 2025 18:02

braghiere added 2 commits October 22, 2025 11:07

Apply JuliaFormatter to fix code style

046bea7

Apply JuliaFormatter to experiment scripts

d239f9a

braghiere requested a review from juliasloan25 October 22, 2025 18:11

braghiere changed the title ~~Rb/spatially varying canopy height~~ spatially varying canopy height Oct 22, 2025

braghiere changed the title ~~spatially varying canopy height~~ Spatially varying canopy height Oct 22, 2025

Fix GPU scalar indexing in canopy height capping tests

6bc905b

Wrap parent() array accesses with Array() to convert from GPU to CPU arrays before scalar indexing. This prevents 'Scalar indexing is disallowed' errors when running tests on GPU.

kmdeck reviewed Oct 22, 2025

View reviewed changes

juliasloan25 reviewed Oct 22, 2025

View reviewed changes

braghiere added 2 commits October 24, 2025 16:29

Refactor spatially-varying canopy height: move capping to preprocessi…

03074f7

…ng, clarify docs, update tests, and clean up experiment scripts

Generalize canopy_height, update tests, and apply JuliaFormatter

839be45

braghiere added 3 commits October 28, 2025 12:13

Remove non-general experiment scripts and Julia tarball for generaliz…

0fc2559

…ation

Move JuliaFormatter to extras/targets to fix Aqua stale_deps failure

10076d5

fix: add JuliaFormatter compat entry for Aqua.jl compliance

25d2c3d

AlexisRenchon reviewed Oct 29, 2025

View reviewed changes

feat: restore backward-compatible API for spatially-varying canopy he…

f4293e1

…ight

		@@ -0,0 +1,187 @@
		# RUN 2 ONLY: Spatially-Varying Canopy Height

Spatially varying canopy height #1502

Are you sure you want to change the base?

Spatially varying canopy height #1502

Conversation

braghiere commented Oct 22, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Purpose

Content

Core Implementation

Testing and Validation

Impact Assessment

To-do

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

kmdeck left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

braghiere commented Oct 29, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

braghiere commented Oct 22, 2025 •

edited

Loading