Marginal Ice advection
This example simulates an (initially) square block of ice moving with a uniform ocean current and deformed by a passing atmospheric cyclone overhead. The scenario mimics a marginal ice zone, where a discrete patch of ice is surrounded by open water, and demonstrates how to:
- set up a marginal ice zone with partial ice cover,
- combine a uniform ocean current with time-varying atmospheric forcing,
- model a translating cyclonic eddy passing over the ice,
- 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:
using ClimaSeaIce
using Oceananigans
using Oceananigans.Units
using Oceananigans.Operators
using Oceananigans.BoundaryConditions
using Printf
using CairoMakie
Lx = 512kilometers
Ly = 256kilometers
Nx = 256
Ny = 128
arch = CPU()
grid = RectilinearGrid(arch; size = (Nx, Ny),
x = (-Lx/2, Lx/2),
y = (0, Ly),
halo = (7, 7),
topology = (Periodic, Bounded, Flat))256×128×1 RectilinearGrid{Float64, Periodic, Bounded, Flat} on CPU with 7×7×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 Boundary conditions
We close the northern and southern walls of the channel with no-slip conditions on the zonal ice velocity:
u_bcs = FieldBoundaryConditions(north = ValueBoundaryCondition(0),
south = ValueBoundaryCondition(0))Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── south: ValueBoundaryCondition: 0
├── north: ValueBoundaryCondition: 0
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)Ocean sea-ice stress
A uniform eastward ocean current provides the background advection that carries the marginal ice block across the channel:
uₒ = 0.2 # m s⁻¹ uniform ocean speed
Uₒ = XFaceField(grid)
set!(Uₒ, uₒ)
fill_halo_regions!(Uₒ)
τₒ = SemiImplicitStress(uₑ=Uₒ)SemiImplicitStress
├── uₑ: 256×128×1 Field{Face, Center, Center} on RectilinearGrid on CPU
├── vₑ: ZeroField{Float64}
├── ρₑ: 1026.0
└── Cᴰ: 0.0055Atmosphere-sea ice stress
The atmosphere consists in a cyclonic eddy translating eastward across the domain. The translation speed is chosen larger than the ocean current so that the cyclone overtakes the ice block during the simulation:
uₐ = 25.0 # m s⁻¹ maximum atmospheric speed modifier
y_cyclone = Ly / 2
@inline x_cyclone(t) = -Lx/2 + 40kilometers + 128kilometers * t / day
@inline cyclone_radius(x, y, t) = sqrt((x - x_cyclone(t))^2 + (y - y_cyclone)^2)
@inline speed(x, y, t) = 1 / 100 * exp(- cyclone_radius(x, y, t) / 80kilometers)
@inline ua_time(x, y, t) = - uₐ * speed(x, y, t) * ( cosd(72) * (x - x_cyclone(t)) + sind(72) * (y - y_cyclone)) / 1000
@inline va_time(x, y, t) = - uₐ * speed(x, y, t) * (- sind(72) * (x - x_cyclone(t)) + cosd(72) * (y - y_cyclone)) / 1000va_time (generic function with 1 method)Initialize the stress at time t = 0:
Uₐ = XFaceField(grid)
Vₐ = YFaceField(grid)
set!(Uₐ, (x, y) -> ua_time(x, y, 0))
set!(Vₐ, (x, y) -> va_time(x, y, 0))
fill_halo_regions!((Uₐ, Vₐ))
τₐu = Field(- Uₐ * sqrt(Uₐ^2 + Vₐ^2) * 1.3 * 1.2e-3)
τₐv = Field(- Vₐ * sqrt(Uₐ^2 + Vₐ^2) * 1.3 * 1.2e-3)256×129×1 Field{Center, Face, Center} on RectilinearGrid on CPU
├── grid: 256×128×1 RectilinearGrid{Float64, Periodic, Bounded, Flat} on CPU with 7×7×0 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: Nothing
├── operand: MultiaryOperation at (Center, Face, Center)
├── status: time=0.0
└── data: 270×143×1 OffsetArray(::Array{Float64, 3}, -6:263, -6:136, 1:1) with eltype Float64 with indices -6:263×-6:136×1:1
└── max=0.0649863, min=-0.0844179, mean=-0.0157795Model configuration
We use an elasto-visco-plastic rheology with split-explicit time stepping and WENO seventh order advection for ice thickness and concentration. We include Coriolis effects appropriate to mid-latitudes:
dynamics = SeaIceMomentumEquation(grid;
top_momentum_stress = (u=τₐu, v=τₐv),
bottom_momentum_stress = τₒ,
coriolis = FPlane(f=1e-4),
rheology = ElastoViscoPlasticRheology(),
solver = SplitExplicitSolver(substeps=150))
model = SeaIceModel(grid;
dynamics,
advection = IncrementalRemapping(),
timestepper = :ForwardEuler,
boundary_conditions = (; u=u_bcs),
ice_thermodynamics = nothing)SeaIceModel{Oceananigans.Architectures.CPU, Oceananigans.Grids.RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 256×128×1 RectilinearGrid{Float64, Periodic, Bounded, Flat} on CPU with 7×7×0 halo
├── timestepper: ForwardEulerTimeStepper
├── ice_thermodynamics: Nothing
├── snow_thermodynamics: Nothing
├── advection: IncrementalRemapping(quadrature_nodes=3)
└── external_heat_fluxes:
├── top: Nothing
└── bottom: Int64Initial conditions
We initialize a square block of ice in the middle of the channel with smoothed edges to limit Gibbs oscillations at the ice edge. The same mask is used for both ice thickness and concentration:
block_radius = 64kilometers
edge_width = 4kilometers
@inline ice_block(x, y) = 0.25 * (1 - tanh((abs(x) - block_radius) / edge_width)) *
(1 - tanh((abs(y - Ly/2) - block_radius) / edge_width))
set!(model, h = ice_block, ℵ = ice_block)Running the simulation
We run the model for 5 days with a 5-minute time step:
simulation = Simulation(model, Δt = 5minutes, stop_time = 5days)Simulation of SeaIceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 5 minutes
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 5 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 entriesTime-varying wind stress
A callback updates the atmospheric stress at each time step as the cyclone translates eastward across the channel:
function compute_wind_stress(sim)
time = sim.model.clock.time
@inline ua(x, y) = ua_time(x, y, time)
@inline va(x, y) = va_time(x, y, time)
set!(Uₐ, ua)
set!(Vₐ, va)
fill_halo_regions!((Uₐ, Vₐ))
compute!(τₐu)
compute!(τₐv)
return nothing
end
simulation.callbacks[:top_stress] = Callback(compute_wind_stress, IterationInterval(1))Callback of compute_wind_stress on IterationInterval(1)Output writer
We save the ice thickness, concentration, and velocity fields:
h = model.ice_thickness
ℵ = model.ice_concentration
u, v = model.velocities
outputs = (; h, u, v, ℵ)
simulation.output_writers[:sea_ice] = JLD2Writer(model, outputs;
filename = "marginal_ice_advection.jld2",
schedule = IterationInterval(10),
overwrite_files = true)JLD2Writer scheduled on IterationInterval(10):
├── filepath: build/literated/marginal_ice_advection.jld2
├── 4 outputs: (h, u, v, ℵ)
├── array_type: Array{Float32}
├── including: [:grid]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)We can finally run the simulation
run!(simulation)[ Info: Initializing simulation...
[ Info: ... simulation initialization complete (40.203 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (5.279 seconds).
[ Info: Simulation is stopping after running for 20.099 minutes.
[ Info: Simulation time 5 days equals or exceeds stop time 5 days.Visualizing the results
We load the saved time series and create an animation:
htimeseries = FieldTimeSeries("marginal_ice_advection.jld2", "h")
utimeseries = FieldTimeSeries("marginal_ice_advection.jld2", "u")
vtimeseries = FieldTimeSeries("marginal_ice_advection.jld2", "v")
ℵtimeseries = FieldTimeSeries("marginal_ice_advection.jld2", "ℵ")
Nt = length(htimeseries)
iter = Observable(1)
𝒱i = @lift(interior(ℵtimeseries[$iter], :, :, 1) .* interior(htimeseries[$iter], :, :, 1))
ℵi = @lift(interior(ℵtimeseries[$iter], :, :, 1))
ui = @lift(interior(utimeseries[$iter], :, :, 1))
vi = @lift(interior(vtimeseries[$iter], :, :, 1))
fig = Figure(size = (1200, 800))
ax1 = Axis(fig[1, 1], title = "Sea ice volume per area, ℵh (m)")
ax2 = Axis(fig[1, 2], title = "Sea ice concentration")
ax3 = Axis(fig[2, 1], title = "Zonal velocity (m s⁻¹)")
ax4 = Axis(fig[2, 2], title = "Meridional velocity (m s⁻¹)")
heatmap!(ax1, 𝒱i, colormap = :magma, colorrange = (0.0, 1.0))
heatmap!(ax2, ℵi, colormap = Reverse(:deep), colorrange = (0.0, 1.0))
heatmap!(ax3, ui, colormap = :balance, colorrange = (-0.2, 0.2))
heatmap!(ax4, vi, colormap = :balance, colorrange = (-0.2, 0.2))
CairoMakie.record(fig, "marginal_ice_advection.mp4", 1:Nt, framerate = 8) do i
iter[] = i
endThe animation shows the ice block being advected eastward by the uniform ocean current while being progressively deformed as the cyclone passes overhead. Ridges and leads form along the ice edges as the cyclone's swirling winds compress and shear the marginal ice cover.
This page was generated using Literate.jl.