diff --git a/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl b/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl index ed93f8c9c..0306aca4b 100644 --- a/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl +++ b/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl @@ -9,6 +9,7 @@ export Radiation, ComponentInterfaces, LatitudeDependentAlbedo, + SeaIceAlbedo, SimilarityTheoryFluxes, MomentumRoughnessLength, ScalarRoughnessLength, @@ -67,6 +68,7 @@ end include("radiation.jl") include("latitude_dependent_albedo.jl") include("tabulated_albedo.jl") +include("sea_ice_albedo.jl") # Turbulent fluxes include("roughness_lengths.jl") diff --git a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl index 1292d6311..5bd715e63 100644 --- a/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl @@ -91,9 +91,18 @@ end Tₛ = convert_to_kelvin(sea_ice_properties.temperature_units, Tₛ) end + # Evaluate state-dependent radiation properties at this grid point. + time = Time(clock.time) + σ = interface_properties.radiation.σ + α = stateindex(interface_properties.radiation.α, i, j, kᴺ, grid, time, CCC) + ϵ = stateindex(interface_properties.radiation.ϵ, i, j, kᴺ, grid, time, CCC) + local_radiation = (; σ, α, ϵ) + local_interface_properties = InterfaceProperties(local_radiation, + interface_properties.specific_humidity_formulation, + interface_properties.temperature_formulation, + interface_properties.velocity_formulation) + # Build thermodynamic and dynamic states in the atmosphere and interface. - # Notation: - # ⋅ 𝒰 ≡ "dynamic" state vector (thermodynamics + reference height + velocity) ℂᵃᵗ = atmosphere_properties.thermodynamics_parameters zᵃᵗ = atmosphere_properties.surface_layer_height # elevation of atmos variables relative to interface @@ -132,7 +141,7 @@ end local_atmosphere_state, local_interior_state, downwelling_radiation, - interface_properties, + local_interface_properties, atmosphere_properties, sea_ice_properties) end diff --git a/src/EarthSystemModels/InterfaceComputations/interface_states.jl b/src/EarthSystemModels/InterfaceComputations/interface_states.jl index 4f1ae6c2d..22b8cb1d2 100644 --- a/src/EarthSystemModels/InterfaceComputations/interface_states.jl +++ b/src/EarthSystemModels/InterfaceComputations/interface_states.jl @@ -16,6 +16,12 @@ struct InterfaceProperties{R, Q, T, V} velocity_formulation :: V end +Adapt.adapt_structure(to, p::InterfaceProperties) = + InterfaceProperties(Adapt.adapt(to, p.radiation), + Adapt.adapt(to, p.specific_humidity_formulation), + Adapt.adapt(to, p.temperature_formulation), + Adapt.adapt(to, p.velocity_formulation)) + ##### ##### Interface specific humidity formulations ##### diff --git a/src/EarthSystemModels/InterfaceComputations/sea_ice_albedo.jl b/src/EarthSystemModels/InterfaceComputations/sea_ice_albedo.jl new file mode 100644 index 000000000..9ffc88417 --- /dev/null +++ b/src/EarthSystemModels/InterfaceComputations/sea_ice_albedo.jl @@ -0,0 +1,133 @@ +""" + SeaIceAlbedo{FT, HI, HS, TS} + +Sea ice albedo parameterization following the CCSM3 scheme (Briegleb et al. 2004). + +Computes broadband albedo as a function of ice thickness, snow depth, and surface +temperature. The scheme blends between bare ice and snow-covered albedos, with +a temperature-dependent reduction near the melting point to implicitly represent +melt pond formation. + +Algorithm: +1. Base cold albedos: bare ice (0.53) and snow-covered (0.82) +2. Temperature reduction within 1C of melting: Δα_ice = 0.075, Δα_snow = 0.10 +3. Thin-ice transition to ocean albedo below h_amin = 0.5 m +4. Snow cover interpolation: full snow albedo at h_snow > h_smin = 0.02 m + +References: +- Briegleb, B.P., C.M. Bitz, E.C. Hunke, W.H. Lipscomb, and M.M. Schramm (2004): + Scientific description of the sea ice component in CCSM3. NCAR Tech Note. +- Briegleb, B.P. and B. Light (2007): NCAR/TN-472+STR. +""" +struct SeaIceAlbedo{FT, HI, HS, TS} + # Cold base albedos (broadband, approx 0.52 * vis + 0.48 * nir) + ice_albedo :: FT # 0.52*0.73 + 0.48*0.33 = 0.538 ≈ 0.54 + snow_albedo :: FT # 0.52*0.96 + 0.48*0.68 = 0.825 ≈ 0.83 + # Melt reduction + ice_melt_reduction :: FT # 0.075 + snow_melt_reduction :: FT # 0.10 + melting_temperature :: FT # 0 C + temperature_range :: FT # 1 C + # Thickness scales + ocean_albedo :: FT # 0.06 + minimum_ice_thickness :: FT # 0.5 m + minimum_snow_depth :: FT # 0.02 m + # References to model fields + ice_thickness :: HI + snow_thickness :: HS + surface_temperature :: TS +end + +Adapt.adapt_structure(to, α::SeaIceAlbedo) = + SeaIceAlbedo(α.ice_albedo, + α.snow_albedo, + α.ice_melt_reduction, + α.snow_melt_reduction, + α.melting_temperature, + α.temperature_range, + α.ocean_albedo, + α.minimum_ice_thickness, + α.minimum_snow_depth, + Adapt.adapt(to, α.ice_thickness), + Adapt.adapt(to, α.snow_thickness), + Adapt.adapt(to, α.surface_temperature)) + +""" + SeaIceAlbedo(ice_thickness, snow_thickness, surface_temperature; + ice_albedo = 0.54, + snow_albedo = 0.83, + ice_melt_reduction = 0.075, + snow_melt_reduction = 0.10, + melting_temperature = 0.0, + temperature_range = 1.0, + ocean_albedo = 0.06, + minimum_ice_thickness = 0.5, + minimum_snow_depth = 0.02) + +Construct a CCSM3 sea ice albedo parameterization. Requires references to the sea ice +model's ice thickness, snow thickness, and surface temperature fields. + +Broadband albedos are approximate averages of the visible and near-IR bands +weighted by solar spectrum (52% visible, 48% near-IR): +- ice_albedo ≈ 0.52 x 0.73 + 0.48 x 0.33 ≈ 0.54 +- snow_albedo ≈ 0.52 x 0.96 + 0.48 x 0.68 ≈ 0.83 +""" +function SeaIceAlbedo(ice_thickness, snow_thickness, surface_temperature; + FT = Float64, + ice_albedo = 0.54, + snow_albedo = 0.83, + ice_melt_reduction = 0.075, + snow_melt_reduction = 0.10, + melting_temperature = 0.0, + temperature_range = 1.0, + ocean_albedo = 0.06, + minimum_ice_thickness = 0.5, + minimum_snow_depth = 0.02) + + return SeaIceAlbedo(convert(FT, ice_albedo), + convert(FT, snow_albedo), + convert(FT, ice_melt_reduction), + convert(FT, snow_melt_reduction), + convert(FT, melting_temperature), + convert(FT, temperature_range), + convert(FT, ocean_albedo), + convert(FT, minimum_ice_thickness), + convert(FT, minimum_snow_depth), + ice_thickness, + snow_thickness, + surface_temperature) +end + +Base.summary(::SeaIceAlbedo{FT}) where FT = "SeaIceAlbedo{$FT}" +Base.show(io::IO, α::SeaIceAlbedo{FT}) where FT = + print(io, "SeaIceAlbedo{$FT}(ice=", α.ice_albedo, + ", snow=", α.snow_albedo, ")") + +@inline function stateindex(α::SeaIceAlbedo, i, j, k, grid, time, loc, args...) + @inbounds hi = α.ice_thickness[i, j, 1] + @inbounds Ts = α.surface_temperature[i, j, 1] + + # Snow thickness: may be nothing (no snow model) + hs = get_snow_thickness(α.snow_thickness, i, j, grid) + + # Temperature-dependent reduction (implicit melt ponds) + Tm = α.melting_temperature + ΔT = α.temperature_range + fT = clamp((Ts - Tm + ΔT) / ΔT, zero(Ts), one(Ts)) + + αi = α.ice_albedo - α.ice_melt_reduction * fT + αs = α.snow_albedo - α.snow_melt_reduction * fT + + # Thin ice → transition to ocean albedo + αo = α.ocean_albedo + fh = clamp(hi / α.minimum_ice_thickness, zero(hi), one(hi)) + αi = αo + (αi - αo) * fh + + # Snow cover blending + fs = clamp(hs / α.minimum_snow_depth, zero(hs), one(hs)) + return fs * αs + (1 - fs) * αi +end + +# Helper to handle nothing snow thickness (no snow model) +@inline get_snow_thickness(hs::Nothing, i, j, grid) = zero(grid) +@inline get_snow_thickness(hs, i, j, grid) = @inbounds hs[i, j, 1] diff --git a/test/test_sea_ice_albedo.jl b/test/test_sea_ice_albedo.jl new file mode 100644 index 000000000..6969cc556 --- /dev/null +++ b/test/test_sea_ice_albedo.jl @@ -0,0 +1,123 @@ +include("runtests_setup.jl") + +using NumericalEarth: stateindex +using NumericalEarth.EarthSystemModels.InterfaceComputations: SeaIceAlbedo, + InterfaceProperties, + Radiation + +using Oceananigans.Units: Time + +@testset "SeaIceAlbedo" begin + for arch in test_architectures + A = typeof(arch) + + grid = RectilinearGrid(arch; + size = (4, 4, 1), + extent = (1, 1, 1), + topology = (Periodic, Periodic, Bounded)) + + hi = Field{Center, Center, Nothing}(grid) + hs = Field{Center, Center, Nothing}(grid) + Ts = Field{Center, Center, Nothing}(grid) + + @testset "Constructor [$A]" begin + α = SeaIceAlbedo(hi, hs, Ts) + @test α isa SeaIceAlbedo{Float64} + @test α.ice_albedo == 0.54 + @test α.snow_albedo == 0.83 + @test α.ocean_albedo == 0.06 + + α32 = SeaIceAlbedo(hi, hs, Ts; FT = Float32) + @test α32 isa SeaIceAlbedo{Float32} + end + + @testset "Nothing snow thickness [$A]" begin + α = SeaIceAlbedo(hi, nothing, Ts) + @test α.snow_thickness === nothing + end + + @testset "Cold thick ice, no snow [$A]" begin + # Thick ice, well below melting: should give cold bare-ice albedo + set!(hi, 2.0) # 2 m (>> minimum_ice_thickness = 0.5) + set!(hs, 0.0) + set!(Ts, -20.0) # well below melting + + α = SeaIceAlbedo(hi, hs, Ts) + time = Time(0.0) + loc = (Center, Center, Center) + + @allowscalar begin + val = stateindex(α, 1, 1, 1, grid, time, loc) + # With Ts = -20, fT = clamp((-20 - 0 + 1)/1, 0, 1) = 0 + # αi = 0.54 - 0 = 0.54, no thin-ice effect (fh=1), no snow (fs=0) + @test val ≈ 0.54 + end + end + + @testset "Cold thick ice, deep snow [$A]" begin + set!(hi, 2.0) + set!(hs, 0.1) # >> minimum_snow_depth = 0.02 + set!(Ts, -20.0) + + α = SeaIceAlbedo(hi, hs, Ts) + time = Time(0.0) + loc = (Center, Center, Center) + + @allowscalar begin + val = stateindex(α, 1, 1, 1, grid, time, loc) + # fs = clamp(0.1/0.02, 0, 1) = 1 → full snow albedo + @test val ≈ 0.83 + end + end + + @testset "Melting ice, no snow [$A]" begin + set!(hi, 2.0) + set!(hs, 0.0) + set!(Ts, 0.0) # at melting point + + α = SeaIceAlbedo(hi, hs, Ts) + time = Time(0.0) + loc = (Center, Center, Center) + + @allowscalar begin + val = stateindex(α, 1, 1, 1, grid, time, loc) + # fT = clamp((0 - 0 + 1)/1, 0, 1) = 1 + # αi = 0.54 - 0.075 = 0.465 + @test val ≈ 0.54 - 0.075 + end + end + + @testset "Thin ice transitions to ocean albedo [$A]" begin + set!(hi, 0.0) # no ice + set!(hs, 0.0) + set!(Ts, -10.0) + + α = SeaIceAlbedo(hi, hs, Ts) + time = Time(0.0) + loc = (Center, Center, Center) + + @allowscalar begin + val = stateindex(α, 1, 1, 1, grid, time, loc) + # fh = clamp(0/0.5, 0, 1) = 0 → pure ocean albedo + @test val ≈ 0.06 + end + end + + @testset "Partial snow cover [$A]" begin + set!(hi, 2.0) + set!(hs, 0.01) # half of minimum_snow_depth + set!(Ts, -20.0) + + α = SeaIceAlbedo(hi, hs, Ts) + time = Time(0.0) + loc = (Center, Center, Center) + + @allowscalar begin + val = stateindex(α, 1, 1, 1, grid, time, loc) + # fs = clamp(0.01/0.02, 0, 1) = 0.5 + # αi = 0.54, αs = 0.83 + @test val ≈ 0.5 * 0.83 + 0.5 * 0.54 + end + end + end +end