diff --git a/experiments/README.md b/experiments/README.md
index 860f131842..c10a3416cf 100644
--- a/experiments/README.md
+++ b/experiments/README.md
@@ -20,7 +20,7 @@ Many of the outputs from the experiments are saved using a non-destructive appro
 output directory of `example_outdir`. When run, if the directory does
 not exist relative to the working directory, then `example_outdir` is created. The results of the latest experiment
 run will be saved into `example_outdir/output_xxx`, where xxxx is an increasing counter. The lower indices store the outputs of previous runs. The latest results are also linked to
-in `example_outdir/output_active`, which can be assumed to always contain the most recent ouput. More details on this style of output directory handeling
+in `example_outdir/output_active`, which can be assumed to always contain the most recent output. More details on this style of output directory handling
 can be found [here](https://clima.github.io/ClimaUtilities.jl/dev/outputpathgenerator/#ActiveLinkStyle-(Non-Destructive))
 
 ### benchmarks
@@ -30,7 +30,7 @@ When a benchmark is run, its outputs are saved using `ActiveLinkStyle` into (ben
 
 ### long_runs
 
-The expirements in experiments/long_runs save its image outputs to (experiment_name)\_longrun\_(device_suffix).
+The experiments in experiments/long_runs save its image outputs to (experiment_name)\_longrun\_(device_suffix).
 All other outputs are saved using `ActiveLinkStyle` to (experiment_name)\_longrun\_(device_suffix)/global_diagnostics.
 
 ### integrated and standalone
diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
index bf3a54805a..cd2807c6dd 100644
--- a/experiments/benchmarks/land.jl
+++ b/experiments/benchmarks/land.jl
@@ -96,8 +96,10 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo,
-        NIR_albedo,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
         f_max,
     ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
@@ -110,9 +112,12 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo = PAR_albedo,
-        NIR_albedo = NIR_albedo,
+        PAR_albedo_dry = PAR_albedo_dry,
+        NIR_albedo_dry = NIR_albedo_dry,
+        PAR_albedo_wet = PAR_albedo_wet,
+        NIR_albedo_wet = NIR_albedo_wet,
     )
+
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s
     runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
@@ -137,6 +142,7 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
 
     # Energy Balance model
     ac_canopy = FT(2.5e3)
+
     # Plant Hydraulics and general plant parameters
     SAI = FT(0.0) # m2/m2
     f_root_to_shoot = FT(3.5)
diff --git a/experiments/integrated/global/global_parameters.jl b/experiments/integrated/global/global_parameters.jl
index 479a56ba0a..1ca486706f 100644
--- a/experiments/integrated/global/global_parameters.jl
+++ b/experiments/integrated/global/global_parameters.jl
@@ -13,8 +13,10 @@ spatially_varying_soil_params =
     K_sat,
     S_s,
     θ_r,
-    PAR_albedo,
-    NIR_albedo,
+    PAR_albedo_dry,
+    NIR_albedo_dry,
+    PAR_albedo_wet,
+    NIR_albedo_wet,
     f_max,
 ) = spatially_varying_soil_params
 soil_params = Soil.EnergyHydrologyParameters(
@@ -27,8 +29,10 @@ soil_params = Soil.EnergyHydrologyParameters(
     K_sat,
     S_s,
     θ_r,
-    PAR_albedo = PAR_albedo,
-    NIR_albedo = NIR_albedo,
+    PAR_albedo_dry = PAR_albedo_dry,
+    NIR_albedo_dry = NIR_albedo_dry,
+    PAR_albedo_wet = PAR_albedo_wet,
+    NIR_albedo_wet = NIR_albedo_wet,
 );
 f_over = FT(3.28) # 1/m
 R_sb = FT(1.484e-4 / 1000) # m/s
diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl
index 0c7fc6b2ef..8fb0daf610 100644
--- a/experiments/long_runs/land.jl
+++ b/experiments/long_runs/land.jl
@@ -97,8 +97,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo,
-        NIR_albedo,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
         f_max,
     ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
@@ -111,9 +113,12 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo = PAR_albedo,
-        NIR_albedo = NIR_albedo,
+        PAR_albedo_dry = PAR_albedo_dry,
+        NIR_albedo_dry = NIR_albedo_dry,
+        PAR_albedo_wet = PAR_albedo_wet,
+        NIR_albedo_wet = NIR_albedo_wet,
     )
+
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s
     runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
diff --git a/experiments/long_runs/land_region.jl b/experiments/long_runs/land_region.jl
index de30e148b3..9be516afc4 100644
--- a/experiments/long_runs/land_region.jl
+++ b/experiments/long_runs/land_region.jl
@@ -100,8 +100,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo,
-        NIR_albedo,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
         f_max,
     ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
@@ -114,8 +116,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo = PAR_albedo,
-        NIR_albedo = NIR_albedo,
+        PAR_albedo_dry = PAR_albedo_dry,
+        NIR_albedo_dry = NIR_albedo_dry,
+        PAR_albedo_wet = PAR_albedo_wet,
+        NIR_albedo_wet = NIR_albedo_wet,
     )
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s
diff --git a/experiments/long_runs/snowy_land.jl b/experiments/long_runs/snowy_land.jl
index 3ccf7dbd54..518a657db4 100644
--- a/experiments/long_runs/snowy_land.jl
+++ b/experiments/long_runs/snowy_land.jl
@@ -99,8 +99,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo,
-        NIR_albedo,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
         f_max,
     ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
@@ -113,8 +115,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo = PAR_albedo,
-        NIR_albedo = NIR_albedo,
+        PAR_albedo_dry = PAR_albedo_dry,
+        NIR_albedo_dry = NIR_albedo_dry,
+        PAR_albedo_wet = PAR_albedo_wet,
+        NIR_albedo_wet = NIR_albedo_wet,
     )
     f_over = FT(3.28) # 1/m
     R_sb = FT(1.484e-4 / 1000) # m/s
diff --git a/experiments/long_runs/soil.jl b/experiments/long_runs/soil.jl
index 85b400a2dd..35c7d30049 100644
--- a/experiments/long_runs/soil.jl
+++ b/experiments/long_runs/soil.jl
@@ -94,8 +94,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo,
-        NIR_albedo,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
         f_max,
     ) = spatially_varying_soil_params
     soil_params = Soil.EnergyHydrologyParameters(
@@ -108,8 +110,10 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         K_sat,
         S_s,
         θ_r,
-        PAR_albedo = PAR_albedo,
-        NIR_albedo = NIR_albedo,
+        PAR_albedo_dry = PAR_albedo_dry,
+        NIR_albedo_dry = NIR_albedo_dry,
+        PAR_albedo_wet = PAR_albedo_wet,
+        NIR_albedo_wet = NIR_albedo_wet,
     )
 
     f_over = FT(3.28) # 1/m
@@ -154,7 +158,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     # The approximation arises because the porosity, residual water fraction,
     # and van Genuchtem parameters are spatially varying but treated constant
     # in solving for equilibrium. Finally, we make a plausible but total guess
-    # for the water table depth using the TOPMODEL parameters. 
+    # for the water table depth using the TOPMODEL parameters.
     function hydrostatic_profile(
         lat::FT,
         z::FT,
diff --git a/ext/CreateParametersExt.jl b/ext/CreateParametersExt.jl
index 1e2c2853ce..ddedcd45ec 100644
--- a/ext/CreateParametersExt.jl
+++ b/ext/CreateParametersExt.jl
@@ -279,7 +279,9 @@ that certain parameters must have the same type (e.g, if a field is
 supplied for porosity, it must be supplied for all other parameters
 defined in the interior of the domain). Some parameters are defined only
 on the surface of the domain (e.g albedo), while other are defined everywhere
-(e.g. porosity). These are indicated with types `F` and `SF`.
+(e.g. porosity). These are indicated with types `F` and `SF`. If both dry/wet albedos
+and general albedos are given as keywords, the dry/wet albedos will override the general
+albedos.
 
 Please see the EnergyHydrologyParameters documentation for a complete list.
 """
@@ -319,12 +321,19 @@ function EnergyHydrologyParameters(
     K_sat::F,
     S_s::F,
     θ_r::F,
-    PAR_albedo::SF = 0.2,
-    NIR_albedo::SF = 0.4,
+    PAR_albedo_dry::SF = nothing,
+    NIR_albedo_dry::SF = nothing,
+    PAR_albedo_wet::SF = nothing,
+    NIR_albedo_wet::SF = nothing,
+    PAR_albedo::SFD = 0.2,
+    NIR_albedo::SFD = 0.4,
+    albedo_calc_top_thickness::TD = 0.07,
     kwargs...,
 ) where {
     F <: Union{<:AbstractFloat, ClimaCore.Fields.Field},
-    SF <: Union{<:AbstractFloat, ClimaCore.Fields.Field},
+    SF <: Union{<:AbstractFloat, ClimaCore.Fields.Field, Nothing},
+    SFD <: Union{<:AbstractFloat, ClimaCore.Fields.Field},
+    TD <: AbstractFloat,
     C,
 }
     earth_param_set = LP.LandParameters(toml_dict)
@@ -386,9 +395,19 @@ function EnergyHydrologyParameters(
     parameters = CP.get_parameter_values(toml_dict, name_map, "Land")
     PSE = typeof(earth_param_set)
     FT = CP.float_type(toml_dict)
-    EnergyHydrologyParameters{FT, F, SF, C, PSE}(;
-        PAR_albedo,
-        NIR_albedo,
+    if isnothing(PAR_albedo_dry)
+        PAR_albedo_dry = FT.(PAR_albedo)
+        NIR_albedo_dry = FT.(NIR_albedo)
+        PAR_albedo_wet = FT.(PAR_albedo)
+        NIR_albedo_wet = FT.(NIR_albedo)
+    end
+    albedo_calc_top_thickness = FT(albedo_calc_top_thickness)
+    EnergyHydrologyParameters{FT, F, typeof(PAR_albedo_dry), C, PSE}(;
+        PAR_albedo_wet,
+        NIR_albedo_wet,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        albedo_calc_top_thickness,
         ν,
         ν_ss_om,
         ν_ss_quartz,
diff --git a/land_longrun_gpu/global_diagnostics/output_active b/land_longrun_gpu/global_diagnostics/output_active
new file mode 120000
index 0000000000..cbcd4c8ca2
--- /dev/null
+++ b/land_longrun_gpu/global_diagnostics/output_active
@@ -0,0 +1 @@
+output_0110
\ No newline at end of file
diff --git a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
index 417cc1a5bf..15df581761 100644
--- a/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
+++ b/lib/ClimaLandSimulations/src/Fluxnet/run_fluxnet.jl
@@ -21,6 +21,12 @@ function run_fluxnet(
     # Now we set up the model. For the soil model, we pick
     # a model type and model args:
     soil_domain = domain.land_domain
+    PAR_albedo =
+        FT(0.5) * params.soil.PAR_albedo_wet +
+        FT(0.5) * params.soil.PAR_albedo_dry
+    NIR_albedo =
+        FT(0.5) * params.soil.NIR_albedo_wet +
+        FT(0.5) * params.soil.NIR_albedo_dry
     soil_ps = Soil.EnergyHydrologyParameters(
         FT;
         ν = params.soil.ν,
@@ -34,8 +40,8 @@ function run_fluxnet(
         z_0m = params.soil.z_0m,
         z_0b = params.soil.z_0b,
         emissivity = params.soil.emissivity,
-        PAR_albedo = params.soil.PAR_albedo,
-        NIR_albedo = params.soil.NIR_albedo,
+        PAR_albedo = PAR_albedo,
+        NIR_albedo = NIR_albedo,
     )
 
     soil_args = (domain = soil_domain, parameters = soil_ps)
diff --git a/src/diagnostics/default_diagnostics.jl b/src/diagnostics/default_diagnostics.jl
index e50187426e..896b755677 100644
--- a/src/diagnostics/default_diagnostics.jl
+++ b/src/diagnostics/default_diagnostics.jl
@@ -177,6 +177,7 @@ function default_diagnostics(
             "lhf",
             "shf",
             "ghf",
+            "salb",
         ]
     elseif output_vars == :short
         soilcanopy_diagnostics = [
diff --git a/src/diagnostics/define_diagnostics.jl b/src/diagnostics/define_diagnostics.jl
index 79f6fa431a..dcc58a07b7 100644
--- a/src/diagnostics/define_diagnostics.jl
+++ b/src/diagnostics/define_diagnostics.jl
@@ -605,6 +605,18 @@ function define_diagnostics!(land_model)
             compute_infiltration!(out, Y, p, t, land_model),
     )
 
+    # Soil albedo
+    add_diagnostic_variable!(
+        short_name = "salb",
+        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.",
+        compute! = (out, Y, p, t) ->
+            compute_soil_albedo!(out, Y, p, t, land_model),
+    )
+
+
     # Soil hydraulic conductivity
     add_diagnostic_variable!(
         short_name = "shc",
diff --git a/src/diagnostics/land_compute_methods.jl b/src/diagnostics/land_compute_methods.jl
index a641c93a30..875f52345d 100644
--- a/src/diagnostics/land_compute_methods.jl
+++ b/src/diagnostics/land_compute_methods.jl
@@ -173,6 +173,20 @@ end
 @diagnostic_compute "soil_net_radiation" Union{SoilCanopyModel, LandModel} p.soil.R_n
 @diagnostic_compute "soil_temperature" Union{SoilCanopyModel, LandModel} p.soil.T
 
+function compute_soil_albedo!(
+    out,
+    Y,
+    p,
+    t,
+    land_model::SoilCanopyModel{FT},
+) where {FT}
+    if isnothing(out)
+        return (p.soil.PAR_albedo .+ p.soil.NIR_albedo) ./ 2
+    else
+        @. out = (p.soil.PAR_albedo + p.soil.NIR_albedo) / 2
+    end
+end
+
 # Soil - Turbulent Fluxes
 @diagnostic_compute "soil_latent_heat_flux" Union{SoilCanopyModel, LandModel} p.soil.turbulent_fluxes.lhf
 @diagnostic_compute "soil_sensible_heat_flux" Union{SoilCanopyModel, LandModel} p.soil.turbulent_fluxes.shf
diff --git a/src/integrated/land.jl b/src/integrated/land.jl
index 39327a7977..24796e7879 100644
--- a/src/integrated/land.jl
+++ b/src/integrated/land.jl
@@ -17,7 +17,7 @@ export LandModel
         snow::SnM
     end
 
-A concrete type of land model used for simulating systems with 
+A concrete type of land model used for simulating systems with
 soil, canopy, snow, soilco2.
 $(DocStringExtensions.FIELDS)
 """
@@ -129,10 +129,8 @@ function LandModel{FT}(;
 
     transpiration = Canopy.PlantHydraulics.DiagnosticTranspiration{FT}()
     ground_conditions =
-        PrognosticGroundConditions{FT, typeof(soil.parameters.PAR_albedo)}(
+        PrognosticGroundConditions{FT, typeof(snow.parameters.α_snow)}(
             snow.parameters.α_snow,
-            soil.parameters.PAR_albedo,
-            soil.parameters.NIR_albedo,
         )
     if :energy in propertynames(canopy_component_args)
         energy_model = canopy_component_types.energy(
@@ -260,7 +258,7 @@ lsm_aux_domain_names(m::LandModel) = (
 A method which makes a function; the returned function
 updates the additional auxiliary variables for the integrated model,
 as well as updates the boundary auxiliary variables for all component
-models. 
+models.
 
 This function is called each ode function evaluation, prior to the tendency function
 evaluation.
@@ -293,7 +291,7 @@ function make_update_boundary_fluxes(
             t,
         )
 
-        # Effective (radiative) land properties 
+        # Effective (radiative) land properties
         set_eff_land_radiation_properties!(
             p,
             land.soil.parameters.earth_param_set,
@@ -396,8 +394,8 @@ function lsm_radiant_energy_fluxes!(
     T_canopy =
         ClimaLand.Canopy.canopy_temperature(canopy.energy, canopy, Y, p, t)
 
-    α_soil_PAR = land.soil.parameters.PAR_albedo
-    α_soil_NIR = land.soil.parameters.NIR_albedo
+    α_soil_PAR = p.soil.PAR_albedo
+    α_soil_NIR = p.soil.NIR_albedo
     ϵ_soil = land.soil.parameters.emissivity
     T_soil = ClimaLand.Domains.top_center_to_surface(p.soil.T)
     α_snow_NIR = land.snow.parameters.α_snow
@@ -579,8 +577,6 @@ struct PrognosticGroundConditions{
     F <: Union{FT, ClimaCore.Fields.Field},
 } <: Canopy.AbstractGroundConditions
     α_snow::FT
-    α_soil_PAR::F
-    α_soil_NIR::F
 end
 
 """
@@ -601,7 +597,7 @@ function Canopy.ground_albedo_PAR(
     p,
     t,
 )
-    return @. (1 - p.snow.snow_cover_fraction) * ground.α_soil_PAR +
+    return @. (1 - p.snow.snow_cover_fraction) * p.soil.PAR_albedo +
               p.snow.snow_cover_fraction * ground.α_snow
 end
 
@@ -623,7 +619,7 @@ function Canopy.ground_albedo_NIR(
     p,
     t,
 )
-    return @. (1 - p.snow.snow_cover_fraction) * ground.α_soil_NIR +
+    return @. (1 - p.snow.snow_cover_fraction) * p.soil.NIR_albedo +
               p.snow.snow_cover_fraction * ground.α_snow
 end
 
@@ -632,7 +628,7 @@ end
 
 Returns the "drivers", or forcing variables, for the LandModel.
 
-These consist of atmospheric and radiative forcing, as well as 
+These consist of atmospheric and radiative forcing, as well as
 soil organic carbon.
 """
 function ClimaLand.get_drivers(model::LandModel)
diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl
index 8f4b04bba0..a4a11f1f4d 100644
--- a/src/integrated/soil_canopy_model.jl
+++ b/src/integrated/soil_canopy_model.jl
@@ -103,11 +103,7 @@ function SoilCanopyModel{FT}(;
     )
 
     transpiration = Canopy.PlantHydraulics.DiagnosticTranspiration{FT}()
-    ground_conditions =
-        PrognosticSoilConditions{typeof(soil.parameters.PAR_albedo)}(
-            soil.parameters.PAR_albedo,
-            soil.parameters.NIR_albedo,
-        )
+    ground_conditions = PrognosticSoilConditions()
     if :energy in propertynames(canopy_component_args)
         energy_model = canopy_component_types.energy(
             canopy_component_args.energy.parameters,
@@ -169,7 +165,7 @@ These include the broadband albedo of the land surface
 `α_sfc`, defined as the ratio of SW_u/SW_d,
 and `T_sfc`, defined as the temperature a blackbody with emissivity
 `ϵ_sfc` would have
-in order to emit the same LW_u as the land surface does. This is called the 
+in order to emit the same LW_u as the land surface does. This is called the
 [effective temperature](https://en.wikipedia.org/wiki/Effective_temperature) in some fields,
 and is not the same as the skin temperature (defined e.g. Equation 7.13 of  Bonan, 2019, Climate Change and Terrestrial Ecosystem Modeling.  DOI: 10.1017/9781107339217).
 """
@@ -252,7 +248,7 @@ function make_update_boundary_fluxes(
             Y,
             t,
         )
-        # Effective (radiative) land properties 
+        # Effective (radiative) land properties
         set_eff_land_radiation_properties!(
             p,
             land.soil.parameters.earth_param_set,
@@ -306,8 +302,8 @@ function lsm_radiant_energy_fluxes!(
     T_canopy =
         ClimaLand.Canopy.canopy_temperature(canopy.energy, canopy, Y, p, t)
 
-    α_soil_PAR = land.soil.parameters.PAR_albedo
-    α_soil_NIR = land.soil.parameters.NIR_albedo
+    α_soil_PAR = p.soil.PAR_albedo
+    α_soil_NIR = p.soil.NIR_albedo
     ϵ_soil = land.soil.parameters.emissivity
     T_soil = ClimaLand.Domains.top_center_to_surface(p.soil.T)
 
@@ -389,26 +385,12 @@ function soil_boundary_fluxes!(
 end
 
 """
-     PrognosticSoilConditions{F <: Union{AbstractFloat, ClimaCore.Fields.Field}} <: Canopy.AbstractGroundConditions
-
-A type of Canopy.AbstractGroundConditions to use when the soil model is prognostic and
-of type `EnergyHydrology`. This is required because the canopy model needs albedo of the ground
-in order to compute its update_aux! function, and that function must only depend on the canopy model.
-
-In the future, we will allocate space for albedo in the cache. In that case, we would *not*
-store them here, twice. `PrognosticSoilConditions` would
-then just be a flag, essentially.
+     PrognosticSoilConditions <: Canopy.AbstractGroundConditions
 
-Note that this struct is linked with the EnergyHydrology model. If we ever had a different
-soil model, we might need to construct a different `PrognosticSoilConditions` because
-the fields may be stored in different places.
+ A type of Canopy.AbstractGroundConditions to use when the soil model is prognostic and
+of type `EnergyHydrology`. `PrognosticSoilConditions` functions as a flag and is used for dispatch
 """
-struct PrognosticSoilConditions{
-    F <: Union{AbstractFloat, ClimaCore.Fields.Field},
-} <: Canopy.AbstractGroundConditions
-    α_PAR::F
-    α_NIR::F
-end
+struct PrognosticSoilConditions <: Canopy.AbstractGroundConditions end
 
 """
     Canopy.ground_albedo_PAR(
@@ -428,7 +410,7 @@ function Canopy.ground_albedo_PAR(
     p,
     t,
 )
-    return ground.α_PAR
+    return p.soil.PAR_albedo
 end
 
 """
@@ -449,10 +431,9 @@ function Canopy.ground_albedo_NIR(
     p,
     t,
 )
-    return ground.α_NIR
+    return p.soil.NIR_albedo
 end
 
-
 function ClimaLand.Soil.sublimation_source(
     prognostic_land_components::Val{(:canopy, :soil, :soilco2)},
     FT,
diff --git a/src/simulations/spatial_parameters.jl b/src/simulations/spatial_parameters.jl
index 35627ff6fe..20c7416c14 100644
--- a/src/simulations/spatial_parameters.jl
+++ b/src/simulations/spatial_parameters.jl
@@ -22,22 +22,22 @@ end
 Reads spatially varying parameters for the canopy, from NetCDF files
 based on CLM and MODIS data, and regrids them to the grid defined by the
 `surface_space` of the Clima simulation. Returns a NamedTuple of ClimaCore
-Fields. 
+Fields.
 
-In particular, this file returns a field for 
+In particular, this file returns a field for
 - clumping index Ω
 - albedo and transmissitivy in PAR and NIR bands
 - leaf angle distribution G function parameter χl
 - Medlyn g1
-- C3 flag 
+- C3 flag
 - VCmax25
 
 The values correspond to the value of the dominant PFT at each point.
 
 The NetCDF files are stored in ClimaArtifacts and more detail on their origin
-is provided there. The keyword arguments `regridder_type` and `extrapolation_bc` 
-affect the regridding by (1) changing how we interpolate to ClimaCore points which 
-are not in the data, and (2) changing how extrapolate to points beyond the range of the 
+is provided there. The keyword arguments `regridder_type` and `extrapolation_bc`
+affect the regridding by (1) changing how we interpolate to ClimaCore points which
+are not in the data, and (2) changing how extrapolate to points beyond the range of the
 data.
 """
 function clm_canopy_parameters(
@@ -166,10 +166,11 @@ end
 Reads spatially varying parameters for the soil model, from NetCDF files
 based on SoilGrids and the van Genuchten data product from Gupta et al 2020,
  and regrids them to the grid defined by the
-`subsurface_space` and `surface_space` of the Clima simulation, as appropriate. 
-Returns a NamedTuple of ClimaCore Fields. 
+`subsurface_space` and `surface_space` of the Clima simulation, as appropriate.
+Returns a NamedTuple of ClimaCore Fields. The albedo parameters are read from
+CLM data.
 
-In particular, this file returns a field for 
+In particular, this file returns a field for
 - (α, n, m) (van Genuchten parameters)
 - Ksat
 - porosity
@@ -180,9 +181,9 @@ In particular, this file returns a field for
 - specific storativity.
 
 The NetCDF files are stored in ClimaArtifacts and more detail on their origin
-is provided there. The keyword arguments `regridder_type` and `extrapolation_bc` 
-affect the regridding by (1) changing how we interpolate to ClimaCore points which 
-are not in the data, and (2) changing how extrapolate to points beyond the range of the 
+is provided there. The keyword arguments `regridder_type` and `extrapolation_bc`
+affect the regridding by (1) changing how we interpolate to ClimaCore points which
+are not in the data, and (2) changing how extrapolate to points beyond the range of the
 data.
 """
 function default_spatially_varying_soil_parameters(
@@ -329,6 +330,25 @@ function default_spatially_varying_soil_parameters(
     # Unsure how to handle two masks
     f_max = oceans_to_value.(f_max, mask, FT(0.0))
     f_max = oceans_to_value.(f_max, soil_params_mask_sfc, FT(0.0))
+    PAR_albedo_dry, NIR_albedo_dry, PAR_albedo_wet, NIR_albedo_wet = map(
+        s -> SpaceVaryingInput(
+            joinpath(
+                ClimaLand.Artifacts.clm_data_folder_path(),
+                "soil_properties_map.nc",
+            ),
+            s,
+            surface_space;
+            regridder_type,
+            regridder_kwargs = (; extrapolation_bc,),
+        ),
+        (
+            "PAR_albedo_dry",
+            "NIR_albedo_dry",
+            "PAR_albedo_wet",
+            "NIR_albedo_wet",
+        ),
+    )
+
     return (;
         ν = ν,
         ν_ss_om = ν_ss_om,
@@ -338,8 +358,10 @@ function default_spatially_varying_soil_parameters(
         K_sat = K_sat,
         S_s = S_s,
         θ_r = θ_r,
-        PAR_albedo = PAR_albedo,
-        NIR_albedo = NIR_albedo,
+        PAR_albedo_wet = PAR_albedo_wet,
+        NIR_albedo_wet = NIR_albedo_wet,
+        PAR_albedo_dry = PAR_albedo_dry,
+        NIR_albedo_dry = NIR_albedo_dry,
         f_max = f_max,
         mask = mask,
     )
diff --git a/src/standalone/Soil/Runoff/Runoff.jl b/src/standalone/Soil/Runoff/Runoff.jl
index 27b48c774b..2afc59b961 100644
--- a/src/standalone/Soil/Runoff/Runoff.jl
+++ b/src/standalone/Soil/Runoff/Runoff.jl
@@ -5,11 +5,7 @@ using ClimaCore.Operators: column_integral_definite!
 using ClimaLand
 import ClimaLand: source!
 using ..ClimaLand.Soil:
-    AbstractSoilSource,
-    AbstractSoilModel,
-    RichardsModel,
-    EnergyHydrology,
-    is_saturated
+    AbstractSoilSource, AbstractSoilModel, RichardsModel, EnergyHydrology
 export TOPMODELRunoff,
     NoRunoff,
     SurfaceRunoff,
@@ -17,7 +13,8 @@ export TOPMODELRunoff,
     TOPMODELSubsurfaceRunoff,
     subsurface_runoff_source,
     topmodel_ss_flux,
-    update_runoff!
+    update_runoff!,
+    is_saturated
 
 
 """
@@ -364,4 +361,15 @@ function topmodel_ss_flux(R_sb::FT, f_over::FT, z∇::FT) where {FT}
     return R_sb * exp(-f_over * z∇)
 end
 
+"""
+    is_saturated(twc::FT, ν::FT) where {FT}
+
+A helper function which can be used to indicate whether a layer of soil is
+saturated based on if the total volumetric water content, `twc` is greater
+than porosity `ν`.
+"""
+function is_saturated(twc::FT, ν::FT) where {FT}
+    return ClimaLand.heaviside(twc - ν)
+end
+
 end
diff --git a/src/standalone/Soil/Soil.jl b/src/standalone/Soil/Soil.jl
index 80fc58bcd2..19f7a9a5f7 100644
--- a/src/standalone/Soil/Soil.jl
+++ b/src/standalone/Soil/Soil.jl
@@ -67,6 +67,7 @@ using Thermodynamics
 using SurfaceFluxes
 using StaticArrays
 import SurfaceFluxes.Parameters as SFP
+import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
 import ClimaLand.Domains: Column, HybridBox, SphericalShell
 import ClimaLand:
     AbstractImExModel,
@@ -176,11 +177,11 @@ append_source(src::Nothing, srcs::Tuple)::Tuple = srcs
 include("./retention_models.jl")
 include("./rre.jl")
 include("./energy_hydrology.jl")
-include("./soil_hydrology_parameterizations.jl")
 include("./soil_heat_parameterizations.jl")
 include("Runoff/Runoff.jl")
 using .Runoff
 include("./boundary_conditions.jl")
+include("./soil_hydrology_parameterizations.jl")
 include("Biogeochemistry/Biogeochemistry.jl")
 using .Biogeochemistry
 end
diff --git a/src/standalone/Soil/boundary_conditions.jl b/src/standalone/Soil/boundary_conditions.jl
index a9a58f7b59..c891dfb4a3 100644
--- a/src/standalone/Soil/boundary_conditions.jl
+++ b/src/standalone/Soil/boundary_conditions.jl
@@ -700,6 +700,10 @@ boundary_vars(bc::AtmosDrivenFluxBC, ::ClimaLand.TopBoundary) = (
     :top_bc,
     :top_bc_wvec,
     :sfc_scratch,
+    :PAR_albedo,
+    :NIR_albedo,
+    :sfc_S_e,
+    :sub_sfc_scratch,
     Runoff.runoff_vars(bc.runoff)...,
 )
 
@@ -717,6 +721,10 @@ boundary_var_domain_names(bc::AtmosDrivenFluxBC, ::ClimaLand.TopBoundary) = (
     :surface,
     :surface,
     :surface,
+    :surface,
+    :surface,
+    :surface,
+    :subsurface,
     Runoff.runoff_var_domain_names(bc.runoff)...,
 )
 """
@@ -742,6 +750,10 @@ boundary_var_types(
     NamedTuple{(:water, :heat), Tuple{FT, FT}},
     ClimaCore.Geometry.WVector{FT},
     FT,
+    FT,
+    FT,
+    FT,
+    FT,
     Runoff.runoff_var_types(bc.runoff, FT)...,
 )
 
diff --git a/src/standalone/Soil/energy_hydrology.jl b/src/standalone/Soil/energy_hydrology.jl
index 74a1bf83e5..d459f80d46 100644
--- a/src/standalone/Soil/energy_hydrology.jl
+++ b/src/standalone/Soil/energy_hydrology.jl
@@ -23,8 +23,8 @@ $(DocStringExtensions.FIELDS)
 """
 Base.@kwdef struct EnergyHydrologyParameters{
     FT <: AbstractFloat,
-    F <: Union{<:AbstractFloat, ClimaCore.Fields.Field},
-    SF <: Union{<:AbstractFloat, ClimaCore.Fields.Field},
+    F <: Union{FT, ClimaCore.Fields.Field},
+    SF <: Union{FT, ClimaCore.Fields.Field},
     C,
     PSE,
 }
@@ -62,10 +62,16 @@ Base.@kwdef struct EnergyHydrologyParameters{
     γ::FT
     "Reference temperature for the viscosity factor"
     γT_ref::FT
-    "Soil PAR Albedo"
-    PAR_albedo::SF
-    "Soil NIR Albedo"
-    NIR_albedo::SF
+    "Soil PAR Albedo dry"
+    PAR_albedo_dry::SF
+    "Soil NIR Albedo dry"
+    NIR_albedo_dry::SF
+    "Soil PAR Albedo wet"
+    PAR_albedo_wet::SF
+    "Soil NIR Albedo wet"
+    NIR_albedo_wet::SF
+    "Thickness of top of soil used in albedo calculations (m)"
+    albedo_calc_top_thickness::FT
     "Soil Emissivity"
     emissivity::FT
     "Roughness length for momentum"
@@ -531,11 +537,22 @@ function ClimaLand.make_update_aux(model::EnergyHydrology)
             κ_sat_unfrozen,
             ρc_ds,
             earth_param_set,
+            PAR_albedo_dry,
+            NIR_albedo_dry,
+            PAR_albedo_wet,
+            NIR_albedo_wet,
         ) = model.parameters
 
         @. p.soil.θ_l =
             volumetric_liquid_fraction(Y.soil.ϑ_l, ν - Y.soil.θ_i, θ_r)
 
+        update_albedo!(
+            model.boundary_conditions.top,
+            p,
+            model.domain,
+            model.parameters,
+        )
+
         @. p.soil.κ = thermal_conductivity(
             model.parameters.κ_dry,
             kersten_number(
@@ -745,9 +762,7 @@ Returns the surface albedo field of the
 `EnergyHydrology` soil model.
 """
 function ClimaLand.surface_albedo(model::EnergyHydrology{FT}, Y, p) where {FT}
-    return @. FT(
-        0.5 * model.parameters.PAR_albedo + 0.5 * model.parameters.NIR_albedo,
-    )
+    return @. (p.soil.PAR_albedo + p.soil.NIR_albedo) / 2
 end
 
 """
@@ -853,7 +868,7 @@ end
 
 Computes turbulent surface fluxes for soil at a point on a surface given
 (1) Surface state conditions (`T_sfc`, `θ_l_sfc`, `θ_i_sfc`)
-(2) Surface properties, such as the topographical height of the surface (`h_sfc`),  
+(2) Surface properties, such as the topographical height of the surface (`h_sfc`),
     the displacement height (`d_sfc`), hydraulic parameters (`hydrology_cm_sfc`,
     `ν_sfc, `θ_r_sfc`, `K_sat_sfc`)
 (4) Atmospheric state conditions (`thermal_state_air`, `u_air`)
@@ -866,14 +881,14 @@ This returns an energy flux and a liquid water volume flux, stored in
 a tuple with self explanatory keys.
 
 Notes:
-In this function, we call SF twice - once for ice, and once for water. 
-We then combine the fluxes from them to get the total SH, LH. 
+In this function, we call SF twice - once for ice, and once for water.
+We then combine the fluxes from them to get the total SH, LH.
 The vapor fluxes are applied to ice and water differently - as either a source term
 or a boundary condition-  so they are kept separate.
 
-For ice, we supply a beta factor and do not need an additional surface resistance, 
-and so we use the output of surface_conditions directly. 
-For water, we do, and it's a little complicated how we handle it. 
+For ice, we supply a beta factor and do not need an additional surface resistance,
+and so we use the output of surface_conditions directly.
+For water, we do, and it's a little complicated how we handle it.
 But once we correct E, we compute SH and LH the same as SurfaceFluxes.jl does.
 """
 function soil_turbulent_fluxes_at_a_point(
@@ -1162,7 +1177,7 @@ end
     ice_fraction(θ_l::FT, θ_i::FT, ν::FT, θ_r::FT)::FT where {FT}
 
 Computes and returns the ice fraction, which is the
-fraction  of the vapor flux that is due to sublimation, and 
+fraction of the vapor flux that is due to sublimation, and
 the fraction of the humidity in the air due to ice, as
 
 f = S_i/(S_i+S_l)
@@ -1205,7 +1220,7 @@ end
                    ) where {FT, EP, C}
 
 Computes the resistance of the top of the soil column to
-water vapor diffusion, as a function of the surface 
+water vapor diffusion, as a function of the surface
 volumetric liquid water fraction `θ_l`, the augmented
 liquid water fraction `ϑ_l`,  the volumetric ice water
 fraction `θ_i`, and other soil parameters.
@@ -1231,7 +1246,7 @@ end
 """
     dry_soil_layer_thickness(S_w::FT, S_c::FT, d_ds::FT)::FT where {FT}
 
-Returns the maximum dry soil layer thickness that can develop under vapor flux; 
+Returns the maximum dry soil layer thickness that can develop under vapor flux;
 this is used when computing the soil resistance to vapor flux according to
 Swenson et al (2012)/Sakaguchi and Zeng (2009).
 """
diff --git a/src/standalone/Soil/soil_hydrology_parameterizations.jl b/src/standalone/Soil/soil_hydrology_parameterizations.jl
index 30e0974bf7..ea1ae5d3a3 100644
--- a/src/standalone/Soil/soil_hydrology_parameterizations.jl
+++ b/src/standalone/Soil/soil_hydrology_parameterizations.jl
@@ -7,7 +7,106 @@ export volumetric_liquid_fraction,
     impedance_factor,
     viscosity_factor,
     dψdϑ,
-    is_saturated
+    update_albedo!
+
+"""
+    albedo_from_moisture(surface_eff_sat::FT, albedo_dry::FT, albedo_wet::FT)
+
+Calculates pointwise albedo for any band as a function of soil effective saturation given
+the dry and wet albedo values for that band.
+"""
+function albedo_from_moisture(
+    surface_eff_sat::FT,
+    albedo_dry::FT,
+    albedo_wet::FT,
+) where {FT}
+    return albedo_dry * (1 - surface_eff_sat) + albedo_wet * surface_eff_sat
+end
+
+
+"""
+    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 larger than the specified `albedo_calc_top_thickness`,
+the water content of the top layer is used in the calclulation. For the PAR and NIR bands,
+
+α_band = α_{band,dry} * (1 - S_e) +  α_{band,wet} * (S_e)
+
+where S_e is the relative soil wetness above some depth, `albedo_calc_top_thickness`. This
+is a modified version of Equation (1) of:
+
+Braghiere, R. K., Wang, Y., Gagné-Landmann, A., Brodrick, P. G., Bloom, A. A., Norton,
+A. J., et al. (2023). The importance of hyperspectral soil albedo information for improving
+Earth system model projections. AGU Advances, 4, e2023AV000910. https://doi.org/10.1029/2023AV000910
+
+where effective saturation is used in place of volumetric soil water content.The dry and wet
+albedo values come from a global soil color map and soil color to albedo map from CLM.
+
+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.
+"""
+function update_albedo!(bc::AtmosDrivenFluxBC, p, soil_domain, model_parameters)
+    (;
+        ν,
+        θ_r,
+        PAR_albedo_dry,
+        NIR_albedo_dry,
+        PAR_albedo_wet,
+        NIR_albedo_wet,
+        albedo_calc_top_thickness,
+    ) = model_parameters
+    S_sfc = p.soil.sfc_S_e
+    FT = eltype(soil_domain.fields.Δz_top)
+    # checks if there is at least 1 layer centered within the top soil depth
+    if minimum(soil_domain.fields.Δz_top) < albedo_calc_top_thickness
+        # ∫H_S_e_dz is the integral of effective saturation from (surface-albedo_calc_top_thickness) to surface
+        ∫H_S_e_dz = p.soil.sfc_S_e
+        # ∫H_dz is integral of 1 from (surface-albedo_calc_top_thickness) to surface
+        ∫H_dz = p.soil.sfc_scratch
+        # zero all centers lower than boundary, set everything above to one
+        @. p.soil.sub_sfc_scratch = ClimaLand.heaviside(
+            albedo_calc_top_thickness + sqrt(eps(FT)),
+            soil_domain.fields.z_sfc - soil_domain.fields.z,
+        )
+        ClimaCore.Operators.column_integral_definite!(
+            ∫H_dz,
+            p.soil.sub_sfc_scratch,
+        )
+        # zeros all effective saturation at levels centered lower than boundary
+        @. p.soil.sub_sfc_scratch =
+            ClimaLand.heaviside(
+                albedo_calc_top_thickness + sqrt(eps(FT)),
+                soil_domain.fields.z_sfc - soil_domain.fields.z,
+            ) * effective_saturation(ν, p.soil.θ_l, θ_r)
+        ClimaCore.Operators.column_integral_definite!(
+            ∫H_S_e_dz,
+            p.soil.sub_sfc_scratch,
+        )
+        @. S_sfc = ∫H_S_e_dz / ∫H_dz
+    else
+        # in the case where no layer is centered above boundary, use the values of the top layer
+        S_sfc .= ClimaLand.Domains.top_center_to_surface(
+            effective_saturation.(ν, p.soil.θ_l, θ_r),
+        )
+    end
+    @. p.soil.PAR_albedo =
+        albedo_from_moisture(S_sfc, PAR_albedo_dry, PAR_albedo_wet)
+    @. p.soil.NIR_albedo =
+        albedo_from_moisture(S_sfc, NIR_albedo_dry, NIR_albedo_wet)
+end
+
+"""
+update_albedo!(bc::AbstractWaterBC, p, soil_domain, model_parameters)
+
+    Does nothing for boundary conditions where albedo is not used
+"""
+function update_albedo!(
+    bc::BC,
+    p,
+    soil_domain,
+    model_parameters,
+) where {BC <: AbstractEnergyHydrologyBC} end
 """
     volumetric_liquid_fraction(ϑ_l::FT, ν_eff::FT, θ_r::FT) where {FT}
 
@@ -67,7 +166,7 @@ end
      approximate_ψ_S_slope(cm::vanGenuchten)
 
 An estimate of the slope of the absolute value of the logψ-logS curve.
-Following Lehmann, Assouline, and Or (2008), we linearize the ψ(S) curve about the inflection point (where d²ψ/dS² = 0, at S = (1+m)^(-m)). 
+Following Lehmann, Assouline, and Or (2008), we linearize the ψ(S) curve about the inflection point (where d²ψ/dS² = 0, at S = (1+m)^(-m)).
 """
 function approximate_ψ_S_slope(cm::vanGenuchten)
     m = cm.m
@@ -87,7 +186,7 @@ end
     ) where {FT}
 
 A point-wise function returning the pressure head in
-variably saturated soil, using the van Genuchten matric potential 
+variably saturated soil, using the van Genuchten matric potential
 if the soil is not saturated, and
 an approximation of the positive pressure in the soil if the
 soil is saturated.
@@ -112,7 +211,7 @@ end
 """
    dψdϑ(cm::vanGenuchten{FT}, ϑ, ν, θ_r, S_s)
 
-Computes and returns the derivative of the pressure head 
+Computes and returns the derivative of the pressure head
 with respect to ϑ for the van Genuchten formulation.
 """
 function dψdϑ(cm::vanGenuchten{FT}, ϑ, ν, θ_r, S_s) where {FT}
@@ -192,7 +291,7 @@ end
 """
    dψdϑ(cm::BrooksCorey{FT}, ϑ, ν, θ_r, S_s)
 
-Computes and returns the derivative of the pressure head 
+Computes and returns the derivative of the pressure head
 with respect to ϑ for the Brooks and Corey formulation.
 """
 function dψdϑ(cm::BrooksCorey{FT}, ϑ, ν, θ_r, S_s) where {FT}
@@ -237,7 +336,7 @@ end
     ) where {FT}
 
 A point-wise function returning the pressure head in
-variably saturated soil, using the Brooks and Corey matric potential 
+variably saturated soil, using the Brooks and Corey matric potential
 if the soil is not saturated, and
 an approximation of the positive pressure in the soil if the
 soil is saturated.
@@ -265,7 +364,7 @@ end
         f_i::FT,
         Ω::FT
     ) where {FT}
-Returns the multiplicative factor reducing conductivity when 
+Returns the multiplicative factor reducing conductivity when
 a fraction of ice `f_i` is present.
 
 Only for use with the `EnergyHydrology` model.
@@ -293,14 +392,3 @@ function viscosity_factor(T::FT, γ::FT, γT_ref::FT) where {FT}
     Theta = FT(exp(factor))
     return Theta
 end
-
-"""
-    is_saturated(twc::FT, ν::FT) where {FT}
-
-A helper function which can be used to indicate whether a layer of soil is 
-saturated based on if the total volumetric water content, `twc` is greater
-than porosity `ν`.
-"""
-function is_saturated(twc::FT, ν::FT) where {FT}
-    return ClimaLand.heaviside(twc - ν)
-end
diff --git a/test/integrated/soil_canopy_lsm.jl b/test/integrated/soil_canopy_lsm.jl
index 22a1516fe1..231d4b2745 100644
--- a/test/integrated/soil_canopy_lsm.jl
+++ b/test/integrated/soil_canopy_lsm.jl
@@ -5,26 +5,202 @@ using ClimaCore
 using ClimaLand
 using ClimaLand.Soil
 using ClimaLand.Canopy
+using Dates
+import ClimaLand.Parameters as LP
+using ClimaLand.Soil.Biogeochemistry
+using ClimaLand.Canopy.PlantHydraulics
+for FT in (Float32, Float64)
+    @testset "PrognosticSoil, FT = $FT" begin
+        # Only care about PAR and NIR albedo values
+        NIR_albedo = FT(0.4)
+        PAR_albedo = FT(0.2)
+        # setup SoilCanopyModel with dummy params
+        earth_param_set = LP.LandParameters(FT)
+        radius = FT(6378.1e3)
+        depth = FT(50)
+        nelements = (50, 10)
+        domain = ClimaLand.Domains.SphericalShell(;
+            radius = radius,
+            depth = depth,
+            nelements = nelements,
+            npolynomial = 1,
+            dz_tuple = FT.((10.0, 0.1)),
+        )
+        hcm = vanGenuchten{FT}(; α = FT(0), n = FT(0))
+        # Radiation
+        start_date = DateTime(2005)
+        atmos, radiation = ClimaLand.prescribed_analytic_forcing(FT)
+        top_bc = ClimaLand.Soil.AtmosDrivenFluxBC(atmos, radiation)
+        zero_water_flux = WaterFluxBC((p, t) -> 0.0)
+        zero_heat_flux = HeatFluxBC((p, t) -> 0.0)
+        boundary_fluxes = (;
+            top = top_bc,
+            bottom = WaterHeatBC(;
+                water = zero_water_flux,
+                heat = zero_heat_flux,
+            ),
+        )
+        soil_params = ClimaLand.Soil.EnergyHydrologyParameters(
+            FT;
+            ν = FT(0.1),
+            ν_ss_om = FT(0.01),
+            ν_ss_quartz = FT(0.01),
+            ν_ss_gravel = FT(0.01),
+            hydrology_cm = hcm,
+            K_sat = FT(0),
+            S_s = FT(0),
+            θ_r = FT(0),
+            NIR_albedo,
+            PAR_albedo,
+            emissivity = FT(0),
+            z_0m = FT(0),
+            z_0b = FT(0),
+        )
+        soil_args = (domain = domain, parameters = soil_params)
+        soilco2_top_bc = Soil.Biogeochemistry.AtmosCO2StateBC()
+        soilco2_bot_bc = Soil.Biogeochemistry.SoilCO2FluxBC((p, t) -> 0.0) # no flux
+        soilco2_sources = ()
+        soilco2_boundary_conditions =
+            (; top = soilco2_top_bc, bottom = soilco2_bot_bc)
+        soilco2_ps = SoilCO2ModelParameters(FT; D_ref = FT(0.0))
+        soilco2_args = (;
+            boundary_conditions = soilco2_boundary_conditions,
+            sources = soilco2_sources,
+            domain = domain,
+            parameters = soilco2_ps,
+        )
+        canopy_component_types = (;
+            autotrophic_respiration = Canopy.AutotrophicRespirationModel{FT},
+            radiative_transfer = Canopy.TwoStreamModel{FT},
+            photosynthesis = Canopy.FarquharModel{FT},
+            conductance = Canopy.MedlynConductanceModel{FT},
+            hydraulics = Canopy.PlantHydraulicsModel{FT},
+            energy = Canopy.BigLeafEnergyModel{FT},
+        )
+        autotrophic_respiration_args =
+            (; parameters = Canopy.AutotrophicRespirationParameters(FT))
+        radiative_transfer_args = (;
+            parameters = Canopy.TwoStreamParameters(
+                FT;
+                Ω = FT(0),
+                α_PAR_leaf = FT(0),
+                τ_PAR_leaf = FT(0),
+                α_NIR_leaf = FT(0),
+                τ_NIR_leaf = FT(0),
+            )
+        )
+        conductance_args =
+            (; parameters = Canopy.MedlynConductanceParameters(FT; g1 = FT(0)))
+        photosynthesis_args = (;
+            parameters = Canopy.FarquharParameters(FT, FT(0); Vcmax25 = FT(0))
+        )
+        ai_parameterization = PlantHydraulics.PrescribedSiteAreaIndex{FT}(
+            TimeVaryingInput(identity),
+            FT(0),
+            FT(0),
+        )
 
-@testset "PrognosticSoil" begin
-    FT = Float64
-    α_PAR = 0.2
-    α_NIR = 0.4
-    soil_driver = ClimaLand.PrognosticSoilConditions{FT}(α_PAR, α_NIR)
-    prognostic_land_components = (:canopy, :soil, :soilco2)
-    @test Canopy.ground_albedo_PAR(
-        Val(prognostic_land_components),
-        soil_driver,
-        nothing,
-        nothing,
-        nothing,
-    ) == 0.2
-    @test Canopy.ground_albedo_NIR(
-        Val(prognostic_land_components),
-        soil_driver,
-        nothing,
-        nothing,
-        nothing,
-    ) == 0.4
+        retention_model = PlantHydraulics.LinearRetentionCurve{FT}(FT(0))
+        conductivity_model = PlantHydraulics.Weibull{FT}(FT(0), FT(0), FT(0))
+        plant_hydraulics_ps = Canopy.PlantHydraulics.PlantHydraulicsParameters(;
+            ai_parameterization = ai_parameterization,
+            ν = FT(0),
+            S_s = FT(0),
+            rooting_depth = FT(0),
+            conductivity_model = conductivity_model,
+            retention_model = retention_model,
+        )
+        plant_hydraulics_args = (
+            parameters = plant_hydraulics_ps,
+            n_stem = 0,
+            n_leaf = 1,
+            compartment_midpoints = [FT(0.0)],
+            compartment_surfaces = [FT(0.0), FT(0.0)],
+        )
+        energy_args = (parameters = Canopy.BigLeafEnergyParameters{FT}(FT(0)),)
+        canopy_component_args = (;
+            autotrophic_respiration = autotrophic_respiration_args,
+            radiative_transfer = radiative_transfer_args,
+            photosynthesis = photosynthesis_args,
+            conductance = conductance_args,
+            hydraulics = plant_hydraulics_args,
+            energy = energy_args,
+        )
+        shared_params =
+            Canopy.SharedCanopyParameters{FT, typeof(earth_param_set)}(
+                FT(0),
+                FT(0),
+                earth_param_set,
+            )
+        canopy_model_args = (;
+            parameters = shared_params,
+            domain = ClimaLand.obtain_surface_domain(domain),
+        )
+        runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(;
+            f_over = FT(0),
+            f_max = fill(FT(0), domain.space.surface),
+            R_sb = FT(0),
+        )
+        Csom = ClimaLand.PrescribedSoilOrganicCarbon{FT}(
+            TimeVaryingInput(identity),
+        )
+        land_input = (
+            atmos = atmos,
+            radiation = radiation,
+            runoff = runoff_model,
+            soil_organic_carbon = Csom,
+        )
+        land_input = (
+            atmos = atmos,
+            radiation = radiation,
+            runoff = runoff_model,
+            soil_organic_carbon = Csom,
+        )
+        model = SoilCanopyModel{FT}(;
+            soilco2_type = Soil.Biogeochemistry.SoilCO2Model{FT},
+            soilco2_args = soilco2_args,
+            land_args = land_input,
+            soil_model_type = Soil.EnergyHydrology{FT},
+            soil_args = soil_args,
+            canopy_component_types = canopy_component_types,
+            canopy_component_args = canopy_component_args,
+            canopy_model_args = canopy_model_args,
+        )
 
+        # initialize model
+        Y, p, coords = initialize(model)
+        # check that albedos have been added to cache
+        @test haskey(p.soil, :PAR_albedo)
+        @test haskey(p.soil, :NIR_albedo)
+        # initialize cache, then check that albedos are set to the correct values
+        set_initial_cache! = make_set_initial_cache(model)
+        set_initial_cache!(p, Y, 0.0)
+        canopy_bc = model.canopy.boundary_conditions
+        @test all(
+            Array(
+                parent(
+                    Canopy.ground_albedo_PAR(
+                        Val(canopy_bc.prognostic_land_components),
+                        canopy_bc.ground,
+                        Y,
+                        p,
+                        0.0,
+                    ),
+                ),
+            ) .== PAR_albedo,
+        )
+        @test all(
+            Array(
+                parent(
+                    Canopy.ground_albedo_NIR(
+                        Val(canopy_bc.prognostic_land_components),
+                        canopy_bc.ground,
+                        Y,
+                        p,
+                        0.0,
+                    ),
+                ),
+            ) .== NIR_albedo,
+        )
+    end
 end
diff --git a/test/integrated/soil_snow.jl b/test/integrated/soil_snow.jl
index 90745dc32c..b579ed756c 100644
--- a/test/integrated/soil_snow.jl
+++ b/test/integrated/soil_snow.jl
@@ -52,8 +52,10 @@ hcm = ClimaLand.Soil.vanGenuchten{FT}(; α = vg_α, n = vg_n)
 ν_ss_quartz = FT(1.0)
 ν_ss_gravel = FT(0.0)
 emissivity = FT(0.99)
-PAR_albedo = FT(0.2)
-NIR_albedo = FT(0.4)
+PAR_albedo_dry = FT(0.2)
+NIR_albedo_dry = FT(0.4)
+PAR_albedo_wet = FT(0.1)
+NIR_albedo_wet = FT(0.2)
 z_0m = FT(0.001)
 z_0b = z_0m
 soil_params = ClimaLand.Soil.EnergyHydrologyParameters(
@@ -66,8 +68,10 @@ soil_params = ClimaLand.Soil.EnergyHydrologyParameters(
     K_sat,
     S_s,
     θ_r,
-    PAR_albedo,
-    NIR_albedo,
+    PAR_albedo_dry,
+    NIR_albedo_dry,
+    PAR_albedo_wet,
+    NIR_albedo_wet,
     emissivity,
     z_0m,
     z_0b,
diff --git a/test/simulations/spatial_parameters.jl b/test/simulations/spatial_parameters.jl
index 8606ae33f3..a5610f8cec 100644
--- a/test/simulations/spatial_parameters.jl
+++ b/test/simulations/spatial_parameters.jl
@@ -40,7 +40,14 @@ spatially_varying_soil_params =
         extrapolation_bc = extrapolation_bc,
     )
 param_names3d = (:ν, :ν_ss_om, :ν_ss_quartz, :ν_ss_gravel, :K_sat, :S_s, :θ_r)
-param_names2d = (:PAR_albedo, :NIR_albedo, :f_max, :mask)
+param_names2d = (
+    :PAR_albedo_dry,
+    :NIR_albedo_dry,
+    :PAR_albedo_wet,
+    :NIR_albedo_wet,
+    :f_max,
+    :mask,
+)
 for p in param_names3d
     @test p ∈ propertynames(spatially_varying_soil_params)
     @test axes(getproperty(spatially_varying_soil_params, p)) ==
diff --git a/test/standalone/Soil/climate_drivers.jl b/test/standalone/Soil/climate_drivers.jl
index 85b2fb7df2..fb9e984e31 100644
--- a/test/standalone/Soil/climate_drivers.jl
+++ b/test/standalone/Soil/climate_drivers.jl
@@ -43,8 +43,6 @@ for FT in (Float32, Float64)
         ν_ss_quartz = FT(1.0)
         ν_ss_gravel = FT(0.0)
         emissivity = FT(0.99)
-        PAR_albedo = FT(0.2)
-        NIR_albedo = FT(0.4)
         z_0m = FT(0.001)
         z_0b = z_0m
         # Radiation
@@ -87,24 +85,30 @@ for FT in (Float32, Float64)
                 heat = zero_heat_flux,
             ),
         )
-        params = ClimaLand.Soil.EnergyHydrologyParameters(
-            FT;
-            ν,
-            ν_ss_om,
-            ν_ss_quartz,
-            ν_ss_gravel,
-            hydrology_cm = hcm,
-            K_sat,
-            S_s,
-            θ_r,
-            PAR_albedo,
-            NIR_albedo,
-            emissivity,
-            z_0m,
-            z_0b,
-        )
 
         for domain in soil_domains
+            NIR_albedo_dry = fill(FT(0.4), domain.space.surface)
+            PAR_albedo_dry = fill(FT(0.2), domain.space.surface)
+            NIR_albedo_wet = fill(FT(0.3), domain.space.surface)
+            PAR_albedo_wet = fill(FT(0.1), domain.space.surface)
+            params = ClimaLand.Soil.EnergyHydrologyParameters(
+                FT;
+                ν,
+                ν_ss_om,
+                ν_ss_quartz,
+                ν_ss_gravel,
+                hydrology_cm = hcm,
+                K_sat,
+                S_s,
+                θ_r,
+                NIR_albedo_dry,
+                PAR_albedo_dry,
+                NIR_albedo_wet,
+                PAR_albedo_wet,
+                emissivity,
+                z_0m,
+                z_0b,
+            )
             model = Soil.EnergyHydrology{FT}(;
                 parameters = params,
                 domain = domain,
@@ -141,6 +145,10 @@ for FT in (Float32, Float64)
                 :top_bc,
                 :top_bc_wvec,
                 :sfc_scratch,
+                :PAR_albedo,
+                :NIR_albedo,
+                :sfc_S_e,
+                :sub_sfc_scratch,
                 :infiltration,
                 :bottom_bc,
                 :bottom_bc_wvec,
@@ -196,8 +204,10 @@ for FT in (Float32, Float64)
             T_sfc = ClimaCore.Fields.zeros(surface_space) .+ FT(280.0)
             @test ClimaLand.surface_emissivity(model, Y, p) == emissivity
             @test ClimaLand.surface_height(model, Y, p) == z_sfc
+            PAR_albedo = p.soil.PAR_albedo
+            NIR_albedo = p.soil.NIR_albedo
             @test ClimaLand.surface_albedo(model, Y, p) ==
-                  PAR_albedo / 2 + NIR_albedo / 2
+                  PAR_albedo ./ 2 .+ NIR_albedo ./ 2
             @test ClimaLand.surface_temperature(model, Y, p, t) == T_sfc
 
             conditions = ClimaLand.turbulent_fluxes(
@@ -256,7 +266,7 @@ for FT in (Float32, Float64)
             #          z_0m,
             #          z_0b,
             #          earth_param_set,
-            #) 
+            #)
         end
     end
 end
diff --git a/test/standalone/Soil/soiltest.jl b/test/standalone/Soil/soiltest.jl
index 5d9f85dcb8..512c0152d6 100644
--- a/test/standalone/Soil/soiltest.jl
+++ b/test/standalone/Soil/soiltest.jl
@@ -232,8 +232,11 @@ for FT in (Float32, Float64)
             FT(7), # Ω
             FT(2.64e-2), # γ
             FT(288), # γT_ref
-            FT(0.2), # NIR_albedo
-            FT(0.2), # PAR_albedo
+            FT(0.2), # NIR_albedo_dry
+            FT(0.2), # PAR_albedo_dry
+            FT(0.2), # NIR_albedo_wet
+            FT(0.2), # PAR_albedo_wet
+            FT(0.035), # top_depth
             FT(0.96), # ϵ
             FT(0.001), # z_0m
             FT(0.01), # z_0b
@@ -429,8 +432,11 @@ for FT in (Float32, Float64)
             FT(7),#Ω
             FT(2.64e-2),#γ
             FT(288),#γT_ref
-            FT(0.2),# NIR_albedo
-            FT(0.2), #PAR_albedo
+            FT(0.2),# NIR_albedo_dry
+            FT(0.2), #PAR_albedo_dry
+            FT(0.2),# NIR_albedo_wet
+            FT(0.2), #PAR_albedo_wet
+            FT(0.035), # top_depth
             FT(0.96),#ϵ
             FT(0.001),# z_0m
             FT(0.01), # z_0b