Better defaults for sea_ice_simulation (#535)

* go for it

* Update src/SeaIceSimulations.jl

* Update SeaIceSimulations.jl

* ice strenght

* some.. salinity

* Update SeaIceSimulations.jl

* Update SeaIceSimulations.jl

* some changes

* ugly workaround

* fix drag

* import

* correction

* bugfix

* another one, probably we should export them?

* update with free drift

* update ClimaSeaIce

Author: simone-silvestri

Project.toml

@@ -44,7 +44,7 @@ ClimaOceanReactantExt = "Reactant"
 Adapt = "4"
 CFTime = "0.1, 0.2"
 CUDA = "4, 5"
-ClimaSeaIce = "0.2.6"
+ClimaSeaIce = "0.3"
 CondaPkg = "0.2.28"
 CubicSplines = "0.2"
 DataDeps = "0.7"

experiments/arctic_simulation.jl

@@ -59,7 +59,7 @@ set!(ocean.model, T=Metadatum(:temperature; dataset),
 ##### A Prognostic Sea-ice model
 #####
 
-using ClimaSeaIce.SeaIceMomentumEquations
+using ClimaSeaIce.SeaIceDynamics
 using ClimaSeaIce.Rheologies
 
 # Remember to pass the SSS as a bottom bc to the sea ice!

src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl

@@ -1,6 +1,6 @@
 using Oceananigans.Operators: Δzᶜᶜᶜ
 using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
-using ClimaSeaIce.SeaIceMomentumEquations: x_momentum_stress, y_momentum_stress
+using ClimaSeaIce.SeaIceDynamics: x_momentum_stress, y_momentum_stress
 
 function compute_sea_ice_ocean_fluxes!(coupled_model)
     ocean = coupled_model.ocean

src/SeaIceSimulations.jl

@@ -17,19 +17,21 @@ using Oceananigans.Operators
 using ClimaSeaIce
 using ClimaSeaIce: SeaIceModel, SlabSeaIceThermodynamics, PhaseTransitions, ConductiveFlux
 using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium
+using ClimaSeaIce.SeaIceDynamics: SplitExplicitSolver, SemiImplicitStress, SeaIceMomentumEquation, StressBalanceFreeDrift
+using ClimaSeaIce.Rheologies: IceStrength, ElastoViscoPlasticRheology
 
 using ClimaOcean.OceanSimulations: Default
 
-function sea_ice_simulation(grid;
+function sea_ice_simulation(grid, ocean=nothing;
                             Δt = 5minutes,
-                            ice_salinity = 0, # psu
+                            ice_salinity = 4, # psu
                             advection = nothing, # for the moment
                             tracers = (),
                             ice_heat_capacity = 2100, # J kg⁻¹ K⁻¹
                             ice_consolidation_thickness = 0.05, # m
                             ice_density = 900, # kg m⁻³
-                            dynamics = nothing,
-                            bottom_heat_boundary_condition = IceWaterThermalEquilibrium(),
+                            dynamics = sea_ice_dynamics(grid, ocean),
+                            bottom_heat_boundary_condition = nothing,
                             phase_transitions = PhaseTransitions(; ice_heat_capacity, ice_density),
                             conductivity = 2, # kg m s⁻³ K⁻¹
                             internal_heat_flux = ConductiveFlux(; conductivity))
@@ -41,6 +43,16 @@ function sea_ice_simulation(grid;
     top_surface_temperature = Field{Center, Center, Nothing}(grid)
     top_heat_boundary_condition = PrescribedTemperature(top_surface_temperature)
 
+    if isnothing(bottom_heat_boundary_condition)
+        if isnothing(ocean)
+            surface_ocean_salinity = 0
+        else
+            kᴺ = size(grid, 3)
+            surface_ocean_salinity = interior(ocean.model.tracers.S, :, :, kᴺ:kᴺ)
+        end
+        bottom_heat_boundary_condition = IceWaterThermalEquilibrium(surface_ocean_salinity)
+    end
+
     ice_thermodynamics = SlabSeaIceThermodynamics(grid;
                                                   internal_heat_flux,
                                                   phase_transitions,
@@ -50,16 +62,12 @@ function sea_ice_simulation(grid;
     bottom_heat_flux = Field{Center, Center, Nothing}(grid)
     top_heat_flux    = Field{Center, Center, Nothing}(grid)
 
-    # top_momentum_stress = (u = Field{Face, Center, Nothing}(grid),
-    #                        v = Field{Center, Face, Nothing}(grid))
-
     # Build the sea ice model
     sea_ice_model = SeaIceModel(grid;
                                 ice_salinity,
                                 advection,
                                 tracers,
                                 ice_consolidation_thickness,
-                                # top_momentum_stress,
                                 ice_thermodynamics,
                                 dynamics,
                                 bottom_heat_flux,
@@ -73,4 +81,40 @@ function sea_ice_simulation(grid;
     return sea_ice
 end
 
-end
+function sea_ice_dynamics(grid, ocean=nothing;
+                          sea_ice_ocean_drag_coefficient = 5.5e-3,
+                          rheology = ElastoViscoPlasticRheology(pressure_formulation = IceStrength()),
+                          coriolis = nothing,
+                          free_drift = nothing,
+                          solver = SplitExplicitSolver(120))
+
+    if isnothing(ocean)
+        SSU = Oceananigans.Fields.ZeroField()
+        SSV = Oceananigans.Fields.ZeroField()
+    else
+        SSU = view(ocean.model.velocities.u, :, :, grid.Nz)
+        SSV = view(ocean.model.velocities.v, :, :, grid.Nz)
+        if isnothing(coriolis)
+            coriolis = ocean.model.coriolis
+        end
+    end
+
+    sea_ice_ocean_drag_coefficient = convert(eltype(grid), sea_ice_ocean_drag_coefficient)
+
+    τo  = SemiImplicitStress(uₑ=SSU, vₑ=SSV, Cᴰ=sea_ice_ocean_drag_coefficient)
+    τua = Field{Face, Center, Nothing}(grid)
+    τva = Field{Center, Face, Nothing}(grid)
+
+    if isnothing(free_drift)
+        free_drift = StressBalanceFreeDrift((u=τua, v=τva), τo)
+    end
+
+    return SeaIceMomentumEquation(grid;
+                                  coriolis,
+                                  top_momentum_stress = (u=τua, v=τva),
+                                  bottom_momentum_stress = τo,
+                                  rheology,
+                                  solver)
+end
+
+end

test/test_ocean_sea_ice_model.jl

@@ -4,7 +4,7 @@ using CUDA
 using Oceananigans.OrthogonalSphericalShellGrids
 using ClimaOcean.OceanSeaIceModels: above_freezing_ocean_temperature!
 using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
-using ClimaSeaIce.SeaIceMomentumEquations
+using ClimaSeaIce.SeaIceDynamics
 using ClimaSeaIce.Rheologies
 
 @inline kernel_melting_temperature(i, j, k, grid, liquidus, S) = @inbounds melting_temperature(liquidus, S[i, j, k])
@@ -81,17 +81,7 @@ using ClimaSeaIce.Rheologies
             ##### Coupled ocean-sea ice and prescribed atmosphere
             #####
 
-            # Adding a sea ice model to the coupled model
-            τua = Field{Face, Center, Nothing}(grid)
-            τva = Field{Center, Face, Nothing}(grid)
-
-            dynamics = SeaIceMomentumEquation(grid;
-                                              coriolis = ocean.model.coriolis,
-                                              top_momentum_stress = (u=τua, v=τva),
-                                              rheology = ElastoViscoPlasticRheology(),
-                                              solver = SplitExplicitSolver(120))
-
-            sea_ice  = sea_ice_simulation(grid; dynamics, advection=WENO(order=7))
+            sea_ice  = sea_ice_simulation(grid, ocean; advection=WENO(order=7))
             liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus
 
             # Set the ocean temperature and salinity

test/test_surface_fluxes.jl

@@ -18,7 +18,7 @@ using Oceananigans.TimeSteppers: update_state!
 using Oceananigans.Units: hours, days
 using ClimaOcean.DataWrangling: all_dates
 
-using ClimaSeaIce.SeaIceMomentumEquations
+using ClimaSeaIce.SeaIceDynamics
 using ClimaSeaIce.Rheologies
 
 import ClimaOcean.OceanSeaIceModels.InterfaceComputations: surface_specific_humidity