Change soil albedo to depend of soil moisture #862

imreddyTeja · 2024-10-17T19:34:38Z

PAR and NIR albedo are moved from the EnergyHydrologyParameters to the auxiliary vars via AtmosDrivenFluxBC. The EnergyHydrologyParameters now hold PAR/NIR wet and dry as a field or float. It also holds top_depth, which is the thickness of the top of the soil used in the albedo calculation.

The constructor for EnergyHydrologyParameters accepts both (PAR/NIR_albedo_dry/wet) and a general PAR/NIR_albedo. Supplying PAR/NIR_albedo is marked as deprecated, and the constructor sets both the dry and wet albedos to the same value.

Previosuly, PAR and NIR albedo were held in a PrognosticSoilConditions
struct, which acted as an intermediary between the soil and canopy models.

Now, PAR and NIR albedo are stored in the cache, and the struct is only
used for dispatch.

Soil albedo is only updated when the top boundary condition is AtmosDrivenFluxBC. Using dispatch for the update_albedo function required moving is_saturated to Runoff.jl and reordering some of the includes in Soil.jl.

The integrated soil-canopy test created a PrognosticSoilConditions
struct and filled in with a PAR and NIR albedo. Then it called
ground_albedo_PAR and ground_albedo_NIR to check the values.

Both of those functions now read from the cache. In order to test them, a
SoilCanopy model would need to be created, initialized, and then updated.
This seems to complex for a test, so this test is deleted.
The long runs and benchmarks now use the wet/dry albedo maps from clm_data. See here

Purpose

Closes #828

To-do

Make a new soil_canopy test
Possibly change calls to EnergyHydrologyParameters constructor

Content

PAR/NIR albedos moved to aux vars via boundry condition
albedo only calculated when AtmosDrivenFluxBC
PAR/NIR albedo depends of top soil moisture volume
EnergyHydrologyParameters holds dry/wet NIR and PAR as float or field
soil albedo added to diagnostics

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

imreddyTeja · 2024-10-22T00:30:57Z

ext/CreateParametersExt.jl

+    ###########DEPRECATED################
+    PAR_albedo::SFD = nothing,
+    NIR_albedo::SFD = nothing,
+    ####################################


Should this be marked as deprecated, or is support going to be continued?

As things are, I think that it would be inconsistent to have both. You set them to be dry/wet, and use the formula to recompute them. So, I think it would end up that the input albedo will be different from the one stored in the cache, which is not desirable.

When the formula recomputes albedo, if both wet and dry are the same, then the calculated value is also that same value. For example if a user passes in PAR_albedo = 0.2, then the constructor sets PAR_albedo_dry = 0.2 and PAR_albedo_wet = 0.2. Then when calculating albedo, the p.soil.PAR_albedo is set to
PAR_albedo_wet * \theta + PAR_albedo_dry * (1-\theta) = 0.2. So isn't the input albedo the same as the one in the cache?

One idea (would have to be a follow-on PR) is to make an AbstractSoilAlbedoModel with two concrete types, one with wet/dry and the other with the bulk (constant) values. Each struct could store their respective values. Then we could also move the albedo parameters outside of EnergyHydrologyParameters and into the AtmosDrivenBC struct.

the top thickness parameter could also live in the struct

the update_albedo function could then dispatch off of BC type and internally on albedo model type.

When the formula recomputes albedo, if both wet and dry are the same, then the calculated value is also that same value. For example if a user passes in PAR_albedo = 0.2, then the constructor sets PAR_albedo_dry = 0.2 and PAR_albedo_wet = 0.2. Then when calculating albedo, the p.soil.PAR_albedo is set to
PAR_albedo_wet * \theta + PAR_albedo_dry * (1-\theta) = 0.2. So isn't the input albedo the same as the one in the cache?

Okay! I thought the integral would change things.

src/standalone/Soil/energy_hydrology.jl

Sbozzolo · 2024-10-22T16:31:51Z

experiments/benchmarks/land.jl

+    PAR_albedo_dry = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "soil_properties_map.nc"),
+        "PAR_albedo_dry",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    NIR_albedo_dry = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "soil_properties_map.nc"),
+        "NIR_albedo_dry",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    PAR_albedo_wet = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "soil_properties_map.nc"),
+        "PAR_albedo_wet",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )
+    NIR_albedo_wet = SpaceVaryingInput(
+        joinpath(clm_artifact_path, "soil_properties_map.nc"),
+        "NIR_albedo_wet",
+        surface_space;
+        regridder_type,
+        regridder_kwargs = (; extrapolation_bc,),
+    )


Can you maybe provide a convenience constructor that sets all of these from a file path?

src/diagnostics/define_diagnostics.jl

src/standalone/Soil/energy_hydrology.jl

Sbozzolo · 2024-10-22T16:36:06Z

ext/CreateParametersExt.jl

+    ###########DEPRECATED################
+    PAR_albedo::SFD = nothing,
+    NIR_albedo::SFD = nothing,
+    ####################################


As things are, I think that it would be inconsistent to have both. You set them to be dry/wet, and use the formula to recompute them. So, I think it would end up that the input albedo will be different from the one stored in the cache, which is not desirable.

src/standalone/Soil/soil_hydrology_parameterizations.jl

src/diagnostics/land_compute_methods.jl

src/standalone/Soil/soil_hydrology_parameterizations.jl

test/integrated/soil_canopy_lsm.jl

kmdeck · 2024-10-22T20:04:15Z

CLM reference: Lawrence, P.J., and Chase, T.N. 2007. Representing a MODIS consistent land surface in the Community Land Model
(CLM 3.0). J. Geophys. Res. 112:G01023. DOI:10.1029/2006JG000168.

src/standalone/Soil/soil_hydrology_parameterizations.jl

kmdeck · 2024-10-22T20:09:55Z

src/standalone/Soil/soil_hydrology_parameterizations.jl

+    ∫H_S_e_dz = p.soil.sfc_θ_l
+    FT = eltype(soil_domain.fields.Δz_top)
+    # checks if there is at least 1 layer centered within the top soil depth
+    if maximum(soil_domain.fields.Δz_top) < FT(0.5) * alebdo_calc_top_thickness


Note that we also have precomputed and stored the layer thickness (face-face):

ClimaLand.jl/src/shared_utilities/Domains.jl

Line 921 in e205018

function get_additional_domain_fields(subsurface_space)

under the name Delta z. You would need to get the top layer thickness with top_center_to_surface.

It doesnt change the computation you do below, but it would allow you to remove the 1/2 everywhere, which might improve readability

hey! I noticed that you changed the code to remove the factor of 1/2, but kept Delta z_top. If we remove the 1/2, we should use Delta z.

Sbozzolo · 2024-10-23T15:14:04Z

src/diagnostics/define_diagnostics.jl

+        long_name = "Soil Albedo",
+        standard_name = "surface albedo",
+        units = "",
+        comments = "Dependent on soil moisture",


Maybe you can add a reference/equation here so that people know how it is calculated? The comment will end up in the NetCDF output

Sbozzolo · 2024-10-23T15:18:32Z

ext/CreateParametersExt.jl

+    ###########DEPRECATED################
+    PAR_albedo::SFD = nothing,
+    NIR_albedo::SFD = nothing,
+    ####################################


When the formula recomputes albedo, if both wet and dry are the same, then the calculated value is also that same value. For example if a user passes in PAR_albedo = 0.2, then the constructor sets PAR_albedo_dry = 0.2 and PAR_albedo_wet = 0.2. Then when calculating albedo, the p.soil.PAR_albedo is set to
PAR_albedo_wet * \theta + PAR_albedo_dry * (1-\theta) = 0.2. So isn't the input albedo the same as the one in the cache?

Okay! I thought the integral would change things.

test/integrated/soil_canopy_lsm.jl

imreddyTeja · 2024-10-24T17:25:12Z

src/diagnostics/define_diagnostics.jl

+        long_name = "Soil Albedo",
+        standard_name = "surface albedo",
+        units = "",
+        comments = "Dependent on soil moisture",


src/standalone/Soil/soil_hydrology_parameterizations.jl

src/diagnostics/land_compute_methods.jl

src/standalone/Soil/soil_hydrology_parameterizations.jl

Sbozzolo · 2024-10-24T17:46:38Z

experiments/benchmarks/land.jl

+    clm_artifact_path = ClimaLand.Artifacts.clm_data_folder_path(; context)
+
+    PAR_albedo_dry, NIR_albedo_dry, PAR_albedo_wet, NIR_albedo_wet =
+        Soil.create_soil_albedo_vars(


I think we can take this one extra step and return direclty the clm data. This function could be clm_soil_albedo_properties and just take the surface_space.

Sbozzolo · 2024-10-24T17:49:32Z

src/diagnostics/land_compute_methods.jl

+function compute_soil_albedo!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel{FT},
+) where {FT}
+    if isnothing(out)
+        return FT(0.5) .* p.soil.PAR_albedo .+ FT(0.5) .* p.soil.NIR_albedo
+    else
+        @. out = FT(0.5) * p.soil.PAR_albedo + FT(0.5) * p.soil.NIR_albedo
+    end
+end


Can this be writen as

@diagnostic_compute "soil_albedo" SoilCanopyModel (p.soil.PAR_albedo .+ p.soil.NIR_albedo) ./ 2

Yes, this works. The docs for diagnostics, however, say
"Or, for convenience, you can use the @diagnostic_compute macro which generates the same function. However, it is better to use that macro only if you are getting a defined variable, such as latent heat flux. (without an operation like the bowen ratio above). "

see docs

Let's stick with you what you have here. I forgot why I added that comment, but I was problably wiser back then.

Sbozzolo · 2024-10-24T17:50:34Z

src/standalone/Soil/energy_hydrology.jl

+    F <: Union{<:FT, ClimaCore.Fields.Field},
+    SF <: Union{<:FT, ClimaCore.Fields.Field},


There's only one subtype of FT ... FT itself, so maybe you want to write F <: Union{FT, ClimaCore.Fields.Field},

Sbozzolo · 2024-10-24T17:51:08Z

src/standalone/Soil/energy_hydrology.jl

-    return @. FT(
-        0.5 * model.parameters.PAR_albedo + 0.5 * model.parameters.NIR_albedo,
-    )
+    return @. FT(0.5 * p.soil.PAR_albedo + 0.5 * p.soil.NIR_albedo)


Suggested change

return @. FT(0.5 * p.soil.PAR_albedo + 0.5 * p.soil.NIR_albedo)

return @. (p.soil.PAR_albedo + p.soil.NIR_albedo) / 2

so that you don't have to worry about the FT

Sbozzolo · 2024-10-24T17:53:08Z

src/standalone/Soil/energy_hydrology.jl

+    PAR_albedo_dry = SpaceVaryingInput(
+        path,
+        "PAR_albedo_dry",
+        surface_space;
+        regridder_type,
+        regridder_kwargs,
+        file_reader_kwargs,
+        compose_function,
+    )
+    NIR_albedo_dry = SpaceVaryingInput(
+        path,
+        "NIR_albedo_dry",
+        surface_space;
+        regridder_type,
+        regridder_kwargs,
+        file_reader_kwargs,
+        compose_function,
+    )
+    PAR_albedo_wet = SpaceVaryingInput(
+        path,
+        "PAR_albedo_wet",
+        surface_space;
+        regridder_type,
+        regridder_kwargs,
+        file_reader_kwargs,
+        compose_function,
+    )
+    NIR_albedo_wet = SpaceVaryingInput(
+        path,
+        "NIR_albedo_wet",
+        surface_space;
+        regridder_type,
+        regridder_kwargs,
+        file_reader_kwargs,
+        compose_function,
+    )
+
+    return PAR_albedo_dry, NIR_albedo_dry, PAR_albedo_wet, NIR_albedo_wet


You could write this as

Suggested change

PAR_albedo_dry = SpaceVaryingInput(

path,

"PAR_albedo_dry",

surface_space;

regridder_type,

regridder_kwargs,

file_reader_kwargs,

compose_function,

)

NIR_albedo_dry = SpaceVaryingInput(

path,

"NIR_albedo_dry",

surface_space;

regridder_type,

regridder_kwargs,

file_reader_kwargs,

compose_function,

)

PAR_albedo_wet = SpaceVaryingInput(

path,

"PAR_albedo_wet",

surface_space;

regridder_type,

regridder_kwargs,

file_reader_kwargs,

compose_function,

)

NIR_albedo_wet = SpaceVaryingInput(

path,

"NIR_albedo_wet",

surface_space;

regridder_type,

regridder_kwargs,

file_reader_kwargs,

compose_function,

)

return PAR_albedo_dry, NIR_albedo_dry, PAR_albedo_wet, NIR_albedo_wet

return map(

s -> SpaceVaryingInput(

path,

s,

surface_space;

regridder_type,

regridder_kwargs,

file_reader_kwargs,

compose_function,

), ("PAR_albedo_dry", "NIR_albedo_dry", "PAR_albedo_wet", "NIT_albedo_dry"))

I also think that we don't need to expose compose_function here. It's hard to see how it would be used.

Sbozzolo · 2024-10-24T17:54:34Z

test/integrated/soil_canopy_lsm.jl

+        atmos = PrescribedAtmosphere(
+            TimeVaryingInput(identity),
+            TimeVaryingInput(identity),
+            TimeVaryingInput(identity),
+            TimeVaryingInput(identity),
+            TimeVaryingInput(identity),
+            TimeVaryingInput(identity),
+            start_date,
+            FT(0),
+            earth_param_set,


You could use prescribed_analytic_forcing to set up radiation and atmosphere

kmdeck · 2024-10-24T20:46:56Z

src/standalone/Soil/soil_hydrology_parameterizations.jl

+the dry and wet albedo values for that band.
+"""
+function albedo_from_moisture(
+    surface_liquid_fraction::FT,


surface_eff_sat::FT

kmdeck · 2024-10-24T20:50:45Z

src/diagnostics/define_diagnostics.jl

+        long_name = "Soil Albedo",
+        standard_name = "surface albedo",
+        units = "",
+        comments = "The mean of PAR and NIR albedo, which are calculated as α_soil_band = α_band_dry * (1 - S_e) + α_band_wet * S_e.",


if you like, you can add this as a default also to the :short diagnostics, and then switch the diagnosticsin the long run back to :short.

kmdeck · 2024-10-24T20:52:02Z

experiments/long_runs/land.jl

@@ -549,7 +555,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
        land,
        start_date;
        output_writer = nc_writer,
-        output_vars = :short,
+        output_vars = :long,


I think you switched to :long in order to get salb saved - but then it is not plotted below (can add it to list of short_names). If you dont want to plot it, we can go back to :short here. Or if you do, we can also add "salb" to the :short diagnostics.

kmdeck · 2024-10-24T21:34:42Z

src/standalone/Soil/soil_hydrology_parameterizations.jl

+    update_albedo!(bc::AtmosDrivenFluxBC, p, soil_domain, model_parameters)
+
+Calculates and updates PAR and NIR albedo as a function of volumetric soil water content at
+the top of the soil. If the soil layers are very large, the water content of the top layer


If the soil layers are larger than a specified parameter alb..., the water content of the top layer...

kmdeck

I think this is ready to merge, but i would like to see the long run output of monthly salb. it would be good to see it has reasonable values first.

kmdeck · 2024-10-24T21:37:35Z

I think this is ready to merge, but i would like to see the long run output of monthly salb. it would be good to see it has reasonable values first.

(this is easy to do with our existing plotting - so if you are unsure - just ask!)

kmdeck · 2024-10-29T15:52:20Z

experiments/benchmarks/bucket.jl

@@ -166,7 +166,7 @@ SciMLBase.solve(prob, ode_algo; dt = Δt, callback = cb)
 @info "Starting profiling"
 # Stop when we profile for MAX_PROFILING_TIME_SECONDS or MAX_PROFILING_SAMPLES
 MAX_PROFILING_TIME_SECONDS = 500
-MAX_PROFILING_SAMPLES = 100
+MAX_PROFILING_SAMPLES = 80


any reason why this change was made?

kmdeck

thank you!

PAR and NIR albedo are moved from the EnergyHydrologyParameters to the auxiliary vars via AtmosDrivenFluxBC. The parameters now hold PAR/NIR wet and dry as a field of float. The old method of creating the parameters with PAR/NIR albedo is marked as deprecated, and the constructor just sets both the wet and dry values to the passed albedo.

update_albedo is no longer called conditionally. Now it is dispatched to a function that either updates albedo or a function that does nothing. update_albedo takes in model.params instead of each param individually. update_albedo now just uses the top layer of if delta_z is larger than the top depth. top_depth is now a parameter instead of hard coded cache var names are changed. Physically irrelevent vars changed to scratch vars Delete integrated soil-canopy test Previosuly, PAR and NIR albedo were held in a PrognosticSoilConditions struct, which acted as an intermediary between the soil and canopy models. Now, PAR and NIR albedo are stored in the cache, and the struct is only used for dispatch. The integrated soil-canopy test created a PrognosticSoilConditions struct and filled in with a PAR and NIR albedo. Then it called ground_albedo_PAR and ground_albedo_NIR to check the values. Both of those functions now read from the cache. In order to test them, a SoilCanopy model would need to be created, initialized, and then updated. This seems to complex for a test, so this test is deleted. Make run_fluxnet work with new SoilEnergy Make update_albedo! comments better Apply suggestions from code review Co-authored-by: Gabriele Bozzola <[email protected]> Change top_depth name to alebdo_calc_top_thickness Add convenience constructor for soil albedo maps Fix calls to create_soil_albedo_vars Fix spelling errors Simplify and correct calculate_albedo! Add citation and fix albedo variable constructor error

This test was deleted because it would not function with previous changes. As requested, this test was recreated to work with the changes. The previous test simply set an albedo it the top bc struct and then checked it using the ground_albedo functions. With the changes, this would only work with an initialized model. Because this is an integrated model test, I thought it made sense to test the SoilCanopy model instead of just creating a Soil model and checking if the function works. This test, like the old version, calls the ground_albedo functions and checks their value. This does seem very involved for a test, so maybe it should just create a EnergyHydrology model and call the function on its Y and p

change create_soil_albedo_vars to clm_soil_albedo_properties and use map in it to reduce code duplication made stylistic change when taking mean of PAR and NIR albedo change <:FT to just FT Move salb to short diagnostics Also make other requested changes such as changing surface_liquid_fraction to surface_eff_sat and making docstring more specific Fix albedo constructor and change benchmark sampling Move salb disgnostic and reduce profiling samples revert samples

Fix merge issue fix merge issue 2

imreddyTeja commented Oct 22, 2024

View reviewed changes

imreddyTeja marked this pull request as ready for review October 22, 2024 16:17

imreddyTeja requested review from kmdeck and Sbozzolo October 22, 2024 16:18

Sbozzolo added Run benchmarks Add this label to run benchmarks on clima Run long runs labels Oct 22, 2024

Sbozzolo reviewed Oct 22, 2024

View reviewed changes