Ice advected on coastline example

This example simulates a solid block of ice moving against a triangular coastline in a periodic channel. The ice is driven by atmospheric winds and interacts with the coastline through immersed boundaries. This example demonstrates how to:

• set up a two-dimensional sea ice model with immersed boundaries,
• 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 two-dimensional periodic channel with a triangular coastline:

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

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

y_max = Ly / 2

arch = CPU()

Oceananigans.Architectures.CPU()

Grid configuration

We create a rectilinear grid with periodic boundaries in x and bounded boundaries in y:

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

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

We define a triangular coastline using an immersed boundary:

bottom(x, y) = ifelse(y > y_max, 0,
               ifelse(abs(x / Lx) * Nx + y / Ly * Ny > 24, 0, 1))

grid = ImmersedBoundaryGrid(grid, GridFittedBoundary(bottom))

256×128×1 ImmersedBoundaryGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded, 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.Bounded, 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}, 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.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Nothing, Nothing, Nothing, @NamedTuple{bottom_and_top::Nothing, south_and_north::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(256, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:256, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_south_and_north_halo!)}, 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.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition}, 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, Bounded, Flat} on CPU with 4×4×0 halo
├── Periodic x ∈ [-256000.0, 256000.0) regularly spaced with Δx=2000.0
├── Bounded  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))

@inline immersed_u_drag(i, j, k, grid, clock, fields, Cᴰ) = @inbounds - Cᴰ * fields.u[i, j, k]
@inline immersed_v_drag(i, j, k, grid, clock, fields, Cᴰ) = @inbounds - Cᴰ * fields.v[i, j, k]

immersed_u_bc = FluxBoundaryCondition(immersed_u_drag, discrete_form=true, parameters=3e-3)
immersed_v_bc = FluxBoundaryCondition(immersed_v_drag, discrete_form=true, parameters=3e-3)

immersed_u_bc = ImmersedBoundaryCondition(top=nothing, bottom=nothing, west=nothing,  east=nothing,  south=immersed_u_bc, north=immersed_u_bc)
immersed_v_bc = ImmersedBoundaryCondition(top=nothing, bottom=nothing, south=nothing, north=nothing, west=immersed_v_bc,  east=immersed_v_bc)

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

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

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: PeriodicBoundaryCondition
├── east: PeriodicBoundaryCondition
├── south: Nothing
├── north: Nothing
├── bottom: Nothing
├── top: Nothing
└── immersed: ImmersedBoundaryCondition with west=Flux, east=Flux, 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 = 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.Bounded, 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
├── 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 (3.315 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (8.718 seconds).
[ Info: Simulation is stopping after running for 25.568 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_on_coastline.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 how the ice block moves and deforms as it interacts with the triangular coastline. The ice accumulates and thickens near the coast due to the no-slip boundary conditions.

This page was generated using Literate.jl.