Implement sea-ice ocean stress (#453)

* add sea-ice ocean stress

* last index

* add this

* add a clock

* add the arctic

* add the metadatum

* thos works!

* bugfix

* Update experiments/arctic_simulation.jl

* Update src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl

* add a test

* remember the allowscalar

* more testing

* fix test

* fix tests

* Update test/test_surface_fluxes.jl

* Update test/test_surface_fluxes.jl

---------

Co-authored-by: Navid C. Constantinou <[email protected]>

Author: simone-silvestri

experiments/arctic_simulation.jl

@@ -52,8 +52,8 @@ ocean = ocean_simulation(grid;
 
 dataset = ECCO4Monthly()
 
-set!(ocean.model, T=Metadata(:temperature; dataset),
-                  S=Metadata(:salinity;    dataset))
+set!(ocean.model, T=Metadatum(:temperature; dataset),
+                  S=Metadatum(:salinity;    dataset))
 
 #####
 ##### A Prognostic Sea-ice model
@@ -63,12 +63,11 @@ using ClimaSeaIce.SeaIceMomentumEquations
 using ClimaSeaIce.Rheologies
 
 # Remember to pass the SSS as a bottom bc to the sea ice!
-SSS = view(ocean.model.tracers.S, :, :, grid.Nz)
+SSS = view(ocean.model.tracers.S.data, :, :, grid.Nz)
 bottom_heat_boundary_condition = IceWaterThermalEquilibrium(SSS)
 
 SSU = view(ocean.model.velocities.u, :, :, grid.Nz)
 SSV = view(ocean.model.velocities.u, :, :, grid.Nz)
-
 τo  = SemiImplicitStress(uₑ=SSU, vₑ=SSV)
 τua = Field{Face, Center, Nothing}(grid)
 τva = Field{Center, Face, Nothing}(grid)
@@ -77,14 +76,14 @@ dynamics = SeaIceMomentumEquation(grid;
                                   coriolis = ocean.model.coriolis,
                                   top_momentum_stress = (u=τua, v=τva),
                                   bottom_momentum_stress = τo,
-                                  ocean_velocities = (u=SSU, v=SSV),
+                                  ocean_velocities = (u=0.1*SSU, v=0.1*SSV),
                                   rheology = ElastoViscoPlasticRheology(),
                                   solver = SplitExplicitSolver(120))
 
 sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics, advection=WENO(order=7))
 
-set!(sea_ice.model, h=Metadata(:sea_ice_thickness;     dataset),
-                    ℵ=Metadata(:sea_ice_concentration; dataset))
+set!(sea_ice.model, h=Metadatum(:sea_ice_thickness;     dataset),
+                    ℵ=Metadatum(:sea_ice_concentration; dataset))
 
 #####
 ##### A Prescribed Atmosphere model
@@ -134,19 +133,20 @@ using Statistics
 
 function progress(sim)
     sea_ice = sim.model.sea_ice
-    hmax = maximum(sea_ice.model.ice_thickness)
-    ℵmax = maximum(sea_ice.model.ice_concentration)
-    hmean = mean(sea_ice.model.ice_thickness)
-    ℵmean = mean(sea_ice.model.ice_concentration)
-    Tmax = maximum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature)
-    Tmin = minimum(sim.model.interfaces.atmosphere_sea_ice_interface.temperature)
+    ocean   = sim.model.ocean
+    hmax  = maximum(sea_ice.model.ice_thickness)
+    ℵmax  = maximum(sea_ice.model.ice_concentration)
+    uimax = maximum(abs, sea_ice.model.velocities.u)
+    vimax = maximum(abs, sea_ice.model.velocities.v)
+    uomax = maximum(abs, ocean.model.velocities.u)
+    vomax = maximum(abs, ocean.model.velocities.v)
 
     step_time = 1e-9 * (time_ns() - wall_time[])
 
     msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt))
     msg2 = @sprintf("max(h): %.2e m, max(ℵ): %.2e ", hmax, ℵmax)
-    msg3 = @sprintf("mean(h): %.2e m, mean(ℵ): %.2e ", hmean, ℵmean)
-    msg4 = @sprintf("extrema(T): (%.2f, %.2f) ᵒC, ", Tmax, Tmin)
+    msg3 = @sprintf("max uᵢ: (%.2f, %.2f) m s⁻¹, ", uimax, vimax)
+    msg4 = @sprintf("max uₒ: (%.2f, %.2f) m s⁻¹, ", uomax, vomax)
     msg5 = @sprintf("wall time: %s \n", prettytime(step_time))
 
     @info msg1 * msg2 * msg3 * msg4 * msg5
@@ -183,16 +183,6 @@ h = FieldTimeSeries("averaged_sea_ice_quantities.jld2", "h")
 hm = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory())
 ℵm = FieldTimeSeries{Center, Center, Nothing}(grid, hE.times; backend=InMemory())
 
-for (i, time) in enumerate(hm.times)
-    counter = 0
-    for (j, t) in enumerate(h.times)
-        if t ≤ time
-            hm[i] .+= h[j]
-            ℵm[i] .+= ℵ[j]
-            counter += 1
-        end
-    end
-    @show counter
-    hm[i] ./= counter
-    ℵm[i] ./= counter
-end
+hMe = [mean(h[t]) for t in 1:length(h.times)]
+ℵMe = [mean(ℵ[t]) for t in 1:length(ℵ.times)]
+

src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl

@@ -4,7 +4,11 @@ using Oceananigans.Operators: ℑxᶠᵃᵃ, ℑyᵃᶠᵃ
 using ClimaOcean.OceanSeaIceModels: sea_ice_concentration
 
 @inline computed_sea_ice_ocean_fluxes(interface) = interface.fluxes
-@inline computed_sea_ice_ocean_fluxes(::Nothing) = (interface_heat = ZeroField(), frazil_heat = ZeroField(), salt = ZeroField())
+@inline computed_sea_ice_ocean_fluxes(::Nothing) = (interface_heat = ZeroField(), 
+                                                    frazil_heat = ZeroField(), 
+                                                    salt = ZeroField(),
+                                                    x_momentum = ZeroField(),
+                                                    y_momentum = ZeroField())
 
 function compute_net_ocean_fluxes!(coupled_model)
     ocean = coupled_model.ocean
@@ -76,8 +80,10 @@ end
     i, j = @index(Global, NTuple)
     kᴺ = size(grid, 3)
     time = Time(clock.time)
-    ρτx = atmos_ocean_fluxes.x_momentum  # zonal momentum flux                      
-    ρτy = atmos_ocean_fluxes.y_momentum  # meridional momentum flux
+    ρτxao = atmos_ocean_fluxes.x_momentum # atmosphere - ocean zonal momentum flux                      
+    ρτyao = atmos_ocean_fluxes.y_momentum # atmosphere - ocean meridional momentum flux
+    ρτxio = sea_ice_ocean_fluxes.x_momentum # sea_ice - ocean zonal momentum flux                      
+    ρτyio = sea_ice_ocean_fluxes.y_momentum # sea_ice - ocean meridional momentum flux
 
     @inbounds begin
         Sₒ = ocean_salinity[i, j, kᴺ]
@@ -121,22 +127,25 @@ end
     ρₒ⁻¹ = 1 / ocean_properties.reference_density
     cₒ   = ocean_properties.heat_capacity
 
-    τxao = ℑxᶠᵃᵃ(i, j, 1, grid, τᶜᶜᶜ, ρₒ⁻¹, ℵ, ρτx)
-    τyao = ℑyᵃᶠᵃ(i, j, 1, grid, τᶜᶜᶜ, ρₒ⁻¹, ℵ, ρτy)
-    Jᵀao = ΣQao  * ρₒ⁻¹ / cₒ
-    Jˢao = - Sₒ * ΣFao
-
-    ρₒ⁻¹ = 1 / ocean_properties.reference_density
-    cₒ   = ocean_properties.heat_capacity
-
     @inbounds begin
         ℵᵢ   = ℵ[i, j, 1]
         Qio  = sea_ice_ocean_fluxes.interface_heat[i, j, 1]
-        Jˢio = sea_ice_ocean_fluxes.salt[i, j, 1]
+
+        Jᵀao = ΣQao  * ρₒ⁻¹ / cₒ
+        Jˢao = - Sₒ * ΣFao
         Jᵀio = Qio * ρₒ⁻¹ / cₒ
+        Jˢio = sea_ice_ocean_fluxes.salt[i, j, 1]
+
+        τxao = ℑxᶠᵃᵃ(i, j, 1, grid, τᶜᶜᶜ, ρₒ⁻¹, ℵ, ρτxao)
+        τyao = ℑyᵃᶠᵃ(i, j, 1, grid, τᶜᶜᶜ, ρₒ⁻¹, ℵ, ρτyao)
+        τxio = ρτxio[i, j, 1] * ρₒ⁻¹ * ℑxᶠᵃᵃ(i, j, 1, grid, ℵ)
+        τyio = ρτyio[i, j, 1] * ρₒ⁻¹ * ℑyᵃᶠᵃ(i, j, 1, grid, ℵ)
+
+        # Stresses
+        τx[i, j, 1] = τxao + τxio
+        τy[i, j, 1] = τyao + τyio
 
-        τx[i, j, 1] = τxao
-        τy[i, j, 1] = τyao
+        # Tracer fluxes
         Jᵀ[i, j, 1] = (1 - ℵᵢ) * Jᵀao + Jᵀio
         Jˢ[i, j, 1] = (1 - ℵᵢ) * Jˢao + Jˢio
     end

src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl

@@ -247,14 +247,18 @@ function sea_ice_ocean_interface(sea_ice::SeaIceSimulation, ocean;
     io_bottom_heat_flux = Field{Center, Center, Nothing}(ocean.model.grid)
     io_frazil_heat_flux = Field{Center, Center, Nothing}(ocean.model.grid)
     io_salt_flux = Field{Center, Center, Nothing}(ocean.model.grid)
+    x_momentum = Field{Face, Center, Nothing}(ocean.model.grid)
+    y_momentum = Field{Center, Face, Nothing}(ocean.model.grid)
 
     @assert io_frazil_heat_flux isa Field{Center, Center, Nothing}
     @assert io_bottom_heat_flux isa Field{Center, Center, Nothing}
     @assert io_salt_flux isa Field{Center, Center, Nothing}
 
     io_fluxes = (interface_heat=io_bottom_heat_flux,
                  frazil_heat=io_frazil_heat_flux,
-                 salt=io_salt_flux)
+                 salt=io_salt_flux,
+                 x_momentum=x_momentum,
+                 y_momentum=y_momentum)
 
     io_properties = (; characteristic_melting_speed)
 

src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl

@@ -1,5 +1,6 @@
 using Oceananigans.Operators: Δzᶜᶜᶜ
 using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
+using ClimaSeaIce.SeaIceMomentumEquations: x_momentum_stress, y_momentum_stress
 
 function compute_sea_ice_ocean_fluxes!(coupled_model)
     ocean = coupled_model.ocean
@@ -18,46 +19,65 @@ function compute_sea_ice_ocean_fluxes!(coupled_model)
 
     ocean_properties = coupled_model.interfaces.ocean_properties
     liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus
-    grid = ocean.model.grid
-    arch = architecture(grid)
+    grid  = ocean.model.grid
+    clock = ocean.model.clock
+    arch  = architecture(grid)
+
+    uᵢ, vᵢ = sea_ice.model.velocities
+    dynamics = sea_ice.model.dynamics
+
+    τs = if isnothing(dynamics)
+        nothing
+    else
+        dynamics.external_momentum_stresses.bottom
+    end
 
     # What about the latent heat removed from the ocean when ice forms?
     # Is it immediately removed from the ocean? Or is it stored in the ice?
-    launch!(arch, grid, :xy, _compute_sea_ice_ocean_latent_heat_flux!,
-            sea_ice_ocean_fluxes, grid, hᵢ, h⁻, ℵᵢ, Sᵢ, Tₒ, Sₒ, liquidus, ocean_properties, interface_properties, Δt)
+    launch!(arch, grid, :xy, _compute_sea_ice_ocean_fluxes!,
+            sea_ice_ocean_fluxes, grid, clock, hᵢ, h⁻, ℵᵢ, Sᵢ, Tₒ, Sₒ, uᵢ, vᵢ, 
+            τs, liquidus, ocean_properties, interface_properties, Δt)
 
     return nothing
 end
 
-@kernel function _compute_sea_ice_ocean_latent_heat_flux!(sea_ice_ocean_fluxes,
-                                                          grid,
-                                                          ice_thickness,
-                                                          previous_ice_thickness,
-                                                          ice_concentration,
-                                                          ice_salinity,
-                                                          ocean_temperature,
-                                                          ocean_salinity,
-                                                          liquidus,
-                                                          ocean_properties,
-                                                          interface_properties,
-                                                          Δt)
-
+@kernel function _compute_sea_ice_ocean_fluxes!(sea_ice_ocean_fluxes,
+                                                grid,
+                                                clock, 
+                                                ice_thickness,
+                                                previous_ice_thickness,
+                                                ice_concentration,
+                                                ice_salinity,
+                                                ocean_temperature,
+                                                ocean_salinity,
+                                                sea_ice_u_velocity, 
+                                                sea_ice_v_velocity, 
+                                                sea_ice_ocean_stresses,
+                                                liquidus,
+                                                ocean_properties,
+                                                interface_properties,
+                                                Δt)
+                                                
     i, j = @index(Global, NTuple)
 
     Nz  = size(grid, 3)
     Qᶠₒ = sea_ice_ocean_fluxes.frazil_heat
     Qᵢₒ = sea_ice_ocean_fluxes.interface_heat
     Jˢ  = sea_ice_ocean_fluxes.salt
+    τx  = sea_ice_ocean_fluxes.x_momentum
+    τy  = sea_ice_ocean_fluxes.y_momentum
+    uᵢ  = sea_ice_u_velocity
+    vᵢ  = sea_ice_v_velocity
     Tₒ  = ocean_temperature
     Sₒ  = ocean_salinity
     Sᵢ  = ice_salinity
     hᵢ  = ice_thickness
+    ℵᵢ  = ice_concentration
     h⁻  = previous_ice_thickness
     ρₒ  = ocean_properties.reference_density
     cₒ  = ocean_properties.heat_capacity
     uₘ★ = interface_properties.characteristic_melting_speed
 
-    ℵ = @inbounds ice_concentration[i, j, 1]
     δQ_frazil = zero(grid)
 
     for k = Nz:-1:1
@@ -94,6 +114,7 @@ end
     @inbounds begin
         Tᴺ = Tₒ[i, j, Nz]
         Sᴺ = Sₒ[i, j, Nz]
+        ℵ  = ℵᵢ[i, j, 1]
     end
 
     # Compute total heat associated with temperature adjustment
@@ -111,6 +132,9 @@ end
     @inbounds Qᶠₒ[i, j, 1] = δQ_frazil
     @inbounds Qᵢₒ[i, j, 1] = δQ_melting * ℵ # Melting depends on concentration
 
+    sea_ice_fields = (; u = uᵢ, v = vᵢ, h = hᵢ, ℵ = ℵᵢ)
+    τₒᵢ = sea_ice_ocean_stresses
+
     @inbounds begin
         # Change in thickness
         Δh = hᵢ[i, j, 1] - h⁻[i, j, 1]
@@ -122,5 +146,9 @@ end
 
         # Update previous ice thickness
         h⁻[i, j, 1] = hᵢ[i, j, 1]
+
+        # Momentum stresses
+        τx[i, j, 1] = x_momentum_stress(i, j, Nz, grid, τₒᵢ, clock, sea_ice_fields)
+        τy[i, j, 1] = y_momentum_stress(i, j, Nz, grid, τₒᵢ, clock, sea_ice_fields)
     end
-end
+end

test/test_surface_fluxes.jl

@@ -186,9 +186,10 @@ end
                                 topology = (Bounded, Bounded, Bounded))
 
         ocean = ocean_simulation(grid; momentum_advection = nothing, 
-                                        tracer_advection = nothing, 
-                                                closure = nothing,
-                                bottom_drag_coefficient = 0.0)
+                                         tracer_advection = nothing, 
+                                                 coriolis = nothing,
+                                                  closure = nothing,
+                                  bottom_drag_coefficient = 0.0)
 
         dates = all_dates(RepeatYearJRA55(), :temperature)
         atmosphere = JRA55PrescribedAtmosphere(arch; end_date=dates[2], backend = InMemory()) 
@@ -210,6 +211,50 @@ end
             # Test that the temperature has snapped up to freezing
             @test minimum(ocean.model.tracers.T) == 0
         end
+
+        @info "Testing Surface Fluxes with sea ice..."
+
+        grid = RectilinearGrid(arch;
+                               size = (2, 2, 2), 
+                             extent = (1, 1, 1),
+                           topology = (Periodic, Periodic, Bounded))
+
+        ocean = ocean_simulation(grid; momentum_advection = nothing, 
+                                         tracer_advection = nothing, 
+                                                 coriolis = nothing,
+                                                  closure = nothing,
+                                  bottom_drag_coefficient = 0.0)
+
+        SSU = view(ocean.model.velocities.u, :, :, grid.Nz)
+        SSV = view(ocean.model.velocities.v, :, :, grid.Nz)
+
+        τo  = SemiImplicitStress(uₑ=SSU, vₑ=SSV, Cᴰ=0.001, ρₑ=1000.0)
+        τua = Field{Face, Center, Nothing}(grid)
+        τva = Field{Center, Face, Nothing}(grid)
+        
+        dynamics = SeaIceMomentumEquation(grid;
+                                          top_momentum_stress = (u=τua, v=τva),
+                                          bottom_momentum_stress = τo,
+                                          rheology = nothing,
+                                          solver = ExplicitSolver())
+
+        sea_ice = sea_ice_simulation(grid; dynamics, advection=Centered())
+
+        # Set a velocity for the ocean
+        fill!(ocean.model.velocities.u, 0.1)
+        fill!(ocean.model.velocities.v, 0.2)
+        fill!(ocean.model.tracers.T,   -2.0)
+
+        # Test that we populate the sea-ice ocean stress
+        earth = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation=Radiation())
+
+        τx = earth.interfaces.sea_ice_ocean_interface.fluxes.x_momentum
+        τy = earth.interfaces.sea_ice_ocean_interface.fluxes.y_momentum
+
+        CUDA.@allowscalar begin
+            @test τx[1, 1, 1] == sqrt(0.1^2 + 0.2^2) * 0.1
+            @test τy[1, 1, 1] == sqrt(0.1^2 + 0.2^2) * 0.2
+        end
     end
 end