Ice advected past an obstacle

This example simulates a solid sheet of ice driven by wind past a rhomboidal island in a doubly periodic domain. The island is represented by an immersed boundary, and the ice sticks to it: the island's faces impose no slip, so the ice piles up on the upwind side and forms a wake downwind. This example demonstrates how to:

• set up a two-dimensional sea ice model with immersed boundaries,
• impose no slip on an immersed boundary with a zero-valued boundary condition,
• prescribe atmospheric and oceanic stresses,
• use elasto-visco-plastic rheology with split-explicit time stepping,
• visualize the evolution of ice thickness, concentration, and velocity.

Install dependencies

First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, ClimaSeaIce, CairoMakie"

The physical domain

We set up a doubly periodic two-dimensional domain holding a single island:

using ClimaSeaIce
using Oceananigans
using Oceananigans.Units
using Oceananigans.Operators
using Printf
using CairoMakie

Lx = 512kilometers
Ly = 256kilometers
Nx = 256
Ny = 128

arch = CPU()

Oceananigans.Architectures.CPU()

Grid configuration

We create a rectilinear grid that is periodic in both $x$ and $y$, so the only boundary the ice can feel is the island itself:

grid = RectilinearGrid(arch; size = (Nx, Ny),
                                x = (-Lx/2, Lx/2),
                                y = (0, Ly),
                             halo = (4, 4),
                         topology = (Periodic, Periodic, Flat))

256×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 4×4×0 halo
├── Periodic x ∈ [-256000.0, 256000.0) regularly spaced with Δx=2000.0
├── Periodic y ∈ [0.0, 256000.0)       regularly spaced with Δy=2000.0
└── Flat z                             

The island is a rhomboid centred in the domain, described by the set of points whose distance from the centre, measured in the $L¹$ norm, is less than one:

island_half_width = 48kilometers
island_half_height = 32kilometers
island_centre = (0, Ly/2)

@inline island(x, y) = abs(x - island_centre[1]) / island_half_width +
                       abs(y - island_centre[2]) / island_half_height < 1

grid = ImmersedBoundaryGrid(grid, GridFittedBoundary(island))

256×128×1 ImmersedBoundaryGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat} on CPU with 4×4×0 halo:
├── immersed_boundary: Oceananigans.ImmersedBoundaries.GridFittedBoundary{Oceananigans.Fields.Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing, Oceananigans.Grids.RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.StaticVerticalDiscretization{Nothing, Nothing, Float64, Float64}, Float64, Float64, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Oceananigans.Architectures.CPU, Oceananigans.Grids.GridSize{256, 128, 1, 4, 4, 0}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Bool, 3, Array{Bool, 3}}, Bool, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, @NamedTuple{bottom_and_top::Nothing, south_and_north::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(256, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:264, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_periodic_south_and_north_halo!)}, 128, 4}, west_and_east::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(136, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:136, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_periodic_west_and_east_halo!)}, 256, 4}}, @NamedTuple{bottom_and_top::Tuple{Nothing, Nothing}, south_and_north::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}, west_and_east::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}}}, Nothing, Nothing}}
├── underlying_grid: 256×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 4×4×0 halo
├── Periodic x ∈ [-256000.0, 256000.0) regularly spaced with Δx=2000.0
├── Periodic y ∈ [0.0, 256000.0)       regularly spaced with Δy=2000.0
└── Flat z                             

Atmospheric and oceanic forcing

We set up atmospheric wind stress. The atmosphere has a constant wind speed:

𝓋ₐ = 10.0   # m s⁻¹
Cᴰ = 1.2e-3 # atmosphere-sea ice drag coefficient
ρₐ = 1.3    # kg m⁻³

1.3

We create a field for the atmospheric wind and compute the stress:

Ua = XFaceField(grid)
τᵤ = Field(- ρₐ * Cᴰ * Ua^2)
set!(Ua, 𝓋ₐ)
compute!(τᵤ)
Oceananigans.BoundaryConditions.fill_halo_regions!(τᵤ)
τᵥ = 0.0

0.0

The ocean stress is represented by a semi-implicit stress with zero ocean velocity:

τₒ = SemiImplicitStress()

SemiImplicitStress
├── uₑ: ZeroField{Float64}
├── vₑ: ZeroField{Float64}
├── ρₑ: 1026.0
└── Cᴰ: 0.0055

Model configuration

We use an elasto-visco-plastic rheology and WENO seventh order for advection of ice thickness and concentration:

dynamics = SeaIceMomentumEquation(grid;
                                  top_momentum_stress = (u=τᵤ, v=τᵥ),
                                  bottom_momentum_stress = τₒ,
                                  rheology = ElastoViscoPlasticRheology(),
                                  solver = SplitExplicitSolver(substeps=150))

SeaIceMomentumEquation
├── coriolis: Nothing
├── rheology: ElastoViscoPlasticRheology{Float64, ClimaSeaIce.Rheologies.ReplacementPressure}
├── auxiliaries: σ₁₁, σ₂₂, σ₁₂, ζᶠᶠᶜ, ζᶜᶜᶜ, Δ, α, uⁿ, vⁿ, P
├── basal_stress: Nothing
├── solver: SplitExplicitSolver{Int64, Symbol}
├── free_drift: nothing
├── external_momentum_stresses: (:top, :bottom)
├── minimum_concentration: 0.001
└── minimum_mass: 1.0

The island holds the ice at rest along its faces. As on a domain boundary, no slip is requested with a zero-valued boundary condition — here in the immersed slot. The rheology reads it and doubles the wall contribution to the shear strain rate, which is the gradient the reflected velocity inside the island would produce. Leaving immersed unset would instead let the ice slide freely along the island.

u_bcs = FieldBoundaryConditions(grid, (Face(), Center(), nothing);
                                immersed = ValueBoundaryCondition(0))

v_bcs = FieldBoundaryConditions(grid, (Center(), Face(), nothing);
                                immersed = ValueBoundaryCondition(0))

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: PeriodicBoundaryCondition
├── east: PeriodicBoundaryCondition
├── south: PeriodicBoundaryCondition
├── north: PeriodicBoundaryCondition
├── bottom: Nothing
├── top: Nothing
└── immersed: ImmersedBoundaryCondition with west=Value, east=Value, south=Nothing, north=Nothing, bottom=Nothing, top=Nothing

We define the model with WENO advection and no thermodynamics:

model = SeaIceModel(grid;
                    advection = WENO(order=7),
                    dynamics,
                    boundary_conditions = (; u=u_bcs, v=v_bcs),
                    ice_thermodynamics = nothing)

SeaIceModel{Oceananigans.Architectures.CPU, Oceananigans.ImmersedBoundaries.ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── grid: 256×128×1 ImmersedBoundaryGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat} on CPU with 4×4×0 halo
├── timestepper: SplitRungeKuttaTimeStepper(3)
├── ice_thermodynamics: Nothing
├── snow_thermodynamics: Nothing
├── advection: WENO{4, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=7)
└── external_heat_fluxes: 
    ├── top: Nothing
    └── bottom: Int64

We initialize the model with uniform ice thickness and concentration:

set!(model, h = 1)
set!(model, ℵ = 1)

Running the simulation

We run the model for 3 days with a 5-minute time step:

simulation = Simulation(model, Δt=5minutes, stop_time=3days)

Simulation of SeaIceModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 5 minutes
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 3 days
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 3 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   └── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
└── output_writers: OrderedDict with no entries

Collecting data

We set up containers to hold the time series data:

htimeseries = []
ℵtimeseries = []
utimeseries = []
vtimeseries = []

function accumulate_timeseries(sim)
    h = sim.model.ice_thickness
    ℵ = sim.model.ice_concentration
    u = sim.model.velocities.u
    v = sim.model.velocities.v
    push!(htimeseries, deepcopy(Array(interior(h))))
    push!(ℵtimeseries, deepcopy(Array(interior(ℵ))))
    push!(utimeseries, deepcopy(Array(interior(u))))
    push!(vtimeseries, deepcopy(Array(interior(v))))
end

simulation.callbacks[:save] = Callback(accumulate_timeseries, IterationInterval(10))

run!(simulation)

[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (531.342 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (9.989 seconds).
[ Info: Simulation is stopping after running for 42.382 minutes.
[ Info: Simulation time 3 days equals or exceeds stop time 3 days.

Visualizing the results

We create an animation of the ice thickness, concentration, and velocity fields:

Nt = length(htimeseries)
iter = Observable(1)

hi = @lift(htimeseries[$iter][:, :, 1])
ℵi = @lift(ℵtimeseries[$iter][:, :, 1])
ui = @lift(utimeseries[$iter][:, :, 1])
vi = @lift(vtimeseries[$iter][:, :, 1])

fig = Figure(size = (1000, 500))
ax = Axis(fig[1, 1], title = "Sea ice thickness (m)")
heatmap!(ax, hi, colormap = :magma, colorrange = (0.0, 2.0))

ax = Axis(fig[1, 2], title = "Sea ice concentration")
heatmap!(ax, ℵi, colormap = Reverse(:deep), colorrange = (0.0, 1))

ax = Axis(fig[2, 1], title = "Zonal velocity (m s⁻¹)")
heatmap!(ax, ui, colorrange = (0, 0.12), colormap = :balance)

ax = Axis(fig[2, 2], title = "Meridional velocity (m s⁻¹)")
heatmap!(ax, vi, colorrange = (-0.025, 0.025), colormap = :bwr)

CairoMakie.record(fig, "sea_ice_advected_past_an_obstacle.mp4", 1:Nt, framerate = 8) do i
    iter[] = i
    @info "Rendering frame $i"
end

[ Info: Rendering frame 1
[ Info: Rendering frame 2
[ Info: Rendering frame 3
[ Info: Rendering frame 4
[ Info: Rendering frame 5
[ Info: Rendering frame 6
[ Info: Rendering frame 7
[ Info: Rendering frame 8
[ Info: Rendering frame 9
[ Info: Rendering frame 10
[ Info: Rendering frame 11
[ Info: Rendering frame 12
[ Info: Rendering frame 13
[ Info: Rendering frame 14
[ Info: Rendering frame 15
[ Info: Rendering frame 16
[ Info: Rendering frame 17
[ Info: Rendering frame 18
[ Info: Rendering frame 19
[ Info: Rendering frame 20
[ Info: Rendering frame 21
[ Info: Rendering frame 22
[ Info: Rendering frame 23
[ Info: Rendering frame 24
[ Info: Rendering frame 25
[ Info: Rendering frame 26
[ Info: Rendering frame 27
[ Info: Rendering frame 28
[ Info: Rendering frame 29
[ Info: Rendering frame 30
[ Info: Rendering frame 31
[ Info: Rendering frame 32
[ Info: Rendering frame 33
[ Info: Rendering frame 34
[ Info: Rendering frame 35
[ Info: Rendering frame 36
[ Info: Rendering frame 37
[ Info: Rendering frame 38
[ Info: Rendering frame 39
[ Info: Rendering frame 40
[ Info: Rendering frame 41
[ Info: Rendering frame 42
[ Info: Rendering frame 43
[ Info: Rendering frame 44
[ Info: Rendering frame 45
[ Info: Rendering frame 46
[ Info: Rendering frame 47
[ Info: Rendering frame 48
[ Info: Rendering frame 49
[ Info: Rendering frame 50
[ Info: Rendering frame 51
[ Info: Rendering frame 52
[ Info: Rendering frame 53
[ Info: Rendering frame 54
[ Info: Rendering frame 55
[ Info: Rendering frame 56
[ Info: Rendering frame 57
[ Info: Rendering frame 58
[ Info: Rendering frame 59
[ Info: Rendering frame 60
[ Info: Rendering frame 61
[ Info: Rendering frame 62
[ Info: Rendering frame 63
[ Info: Rendering frame 64
[ Info: Rendering frame 65
[ Info: Rendering frame 66
[ Info: Rendering frame 67
[ Info: Rendering frame 68
[ Info: Rendering frame 69
[ Info: Rendering frame 70
[ Info: Rendering frame 71
[ Info: Rendering frame 72
[ Info: Rendering frame 73
[ Info: Rendering frame 74
[ Info: Rendering frame 75
[ Info: Rendering frame 76
[ Info: Rendering frame 77
[ Info: Rendering frame 78
[ Info: Rendering frame 79
[ Info: Rendering frame 80
[ Info: Rendering frame 81
[ Info: Rendering frame 82
[ Info: Rendering frame 83
[ Info: Rendering frame 84
[ Info: Rendering frame 85
[ Info: Rendering frame 86
[ Info: Rendering frame 87

The animation shows the ice sheet deforming as it is driven past the island. Ice converges and thickens on the upwind face, and the no-slip condition drags it to a halt along the island's flanks, leaving a slow wake that persists downwind.

This page was generated using Literate.jl.