Sea ice advected by an atmospheric anticyclone example

This example simulates sea ice advected by an atmospheric anticyclone, based on the experiment found in the paper:

Simulating Linear Kinematic Features in Viscous-Plastic Sea Ice Models on Quadrilateral and Triangular Grids With Different Variable Staggering

The ice is driven by a moving anticyclonic atmospheric eddy and interacts with a cyclonic ocean eddy. This example demonstrates how to:

• set up a two-dimensional sea ice model with time-varying atmospheric forcing,
• prescribe moving atmospheric and ocean eddies,
• use elasto-visco-plastic rheology with Coriolis effects,
• save and visualize time series data.

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 bounded domain:

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

arch = CPU()

L  = 512kilometers
𝓋ₒ = 0.01 # m s⁻¹ maximum ocean speed
𝓋ₐ = 30.0 # m s⁻¹ maximum atmospheric speed modifier

grid = RectilinearGrid(arch;
                       size = (128, 128),
                          x = (0, L),
                          y = (0, L),
                       halo = (7, 7),
                   topology = (Bounded, Bounded, Flat))

128×128×1 RectilinearGrid{Float64, Bounded, Bounded, Flat} on CPU with 7×7×0 halo
├── Bounded  x ∈ [0.0, 512000.0] regularly spaced with Δx=4000.0
├── Bounded  y ∈ [0.0, 512000.0] regularly spaced with Δy=4000.0
└── Flat z                       

Boundary conditions

We set value boundary conditions for velocities:

u_bcs = FieldBoundaryConditions(north=ValueBoundaryCondition(0),
                                south=ValueBoundaryCondition(0))

v_bcs = FieldBoundaryConditions(west=ValueBoundaryCondition(0),
                                east=ValueBoundaryCondition(0))

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: ValueBoundaryCondition: 0
├── east: ValueBoundaryCondition: 0
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

Ocean sea-ice stress

We set up constant ocean velocities corresponding to a cyclonic eddy:

Uₒ = XFaceField(grid)
Vₒ = YFaceField(grid)

set!(Uₒ, (x, y) -> 𝓋ₒ * (2y - L) / L)
set!(Vₒ, (x, y) -> 𝓋ₒ * (L - 2x) / L)
fill_halo_regions!((Uₒ, Vₒ))

τₒ = SemiImplicitStress(uₑ=Uₒ, vₑ=Vₒ)

SemiImplicitStress
├── uₑ: 129×128×1 Field{Oceananigans.Grids.Face, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CPU
├── vₑ: 128×129×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CPU
├── ρₑ: 1026.0
└── Cᴰ: 0.0055

Atmosphere-sea ice stress

We set up atmospheric velocities corresponding to an anticyclonic eddy moving northeast. The eddy center moves with time:

@inline center(t) = 256kilometers + 51.2kilometers * t / day
@inline radius(x, y, t) = sqrt((x - center(t))^2 + (y - center(t))^2)
@inline speed(x, y, t) = 1 / 100 * exp(- radius(x, y, t) / 100kilometers)

@inline ua_time(x, y, t) = - 𝓋ₐ * speed(x, y, t) * (  cosd(72) * (x - center(t)) + sind(72) * (y - center(t))) / 1000
@inline va_time(x, y, t) = - 𝓋ₐ * speed(x, y, t) * (- sind(72) * (x - center(t)) + cosd(72) * (y - center(t))) / 1000

va_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)

128×129×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CPU
├── grid: 128×128×1 RectilinearGrid{Float64, Bounded, Bounded, Flat} on CPU with 7×7×0 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: ZeroFlux, east: ZeroFlux, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: Nothing
├── operand: MultiaryOperation at (Center, Face, Center)
├── status: time=0.0
└── data: 142×143×1 OffsetArray(::Array{Float64, 3}, -6:135, -6:136, 1:1) with eltype Float64 with indices -6:135×-6:136×1:1
    └── max=0.189868, min=-0.189868, mean=2.37482e-17

Model configuration

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

dynamics = SeaIceMomentumEquation(grid;
                                  top_momentum_stress = (u=τₐu, v=τₐv),
                                  bottom_momentum_stress = τₒ,
                                  coriolis = FPlane(f=1e-4))

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

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

Initial conditions

We start with a concentration of ℵ = 1 and an initial height field with perturbations around 0.3 meters:

h₀(x, y) = 0.3 + 0.005 * (sin(60 * x / 1000kilometers) + sin(30 * y / 1000kilometers))

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

Running the simulation

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

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

Simulation of SeaIceModel
├── Next time step: 2 minutes
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 2 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

Time-varying wind stress

We need to evolve the wind stress field in time. We set up a callback to update the atmospheric stress at each time step:

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 set up an output writer to 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 = "sea_ice_advected_by_anticyclone.jld2",
                                                 schedule = IterationInterval(5),
                                                 overwrite_existing = true)

JLD2Writer scheduled on IterationInterval(5):
├── filepath: build/literated/sea_ice_advected_by_anticyclone.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 (34.742 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (8.104 seconds).
[ Info: Simulation is stopping after running for 24.216 minutes.
[ Info: Simulation time 2 days equals or exceeds stop time 2 days.

Visualizing the results

We load the saved time series and create an animation:

htimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "h")
utimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "u")
vtimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "v")
ℵtimeseries = FieldTimeSeries("sea_ice_advected_by_anticyclone.jld2", "ℵ")

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

hi = @lift(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, 1200))
ax1 = Axis(fig[1, 1], title = "Sea ice thickness (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, hi, colormap = :magma, colorrange = (0.23, 0.37))
heatmap!(ax2, ℵi, colormap = Reverse(:deep), colorrange = (0.9, 1))
heatmap!(ax3, ui, colorrange = (-0.1, 0.1), colormap = :balance)
heatmap!(ax4, vi, colorrange = (-0.1, 0.1), colormap = :balance)

CairoMakie.record(fig, "sea_ice_advected_by_anticyclone.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
[ Info: Rendering frame 88
[ Info: Rendering frame 89
[ Info: Rendering frame 90
[ Info: Rendering frame 91
[ Info: Rendering frame 92
[ Info: Rendering frame 93
[ Info: Rendering frame 94
[ Info: Rendering frame 95
[ Info: Rendering frame 96
[ Info: Rendering frame 97
[ Info: Rendering frame 98
[ Info: Rendering frame 99
[ Info: Rendering frame 100
[ Info: Rendering frame 101
[ Info: Rendering frame 102
[ Info: Rendering frame 103
[ Info: Rendering frame 104
[ Info: Rendering frame 105
[ Info: Rendering frame 106
[ Info: Rendering frame 107
[ Info: Rendering frame 108
[ Info: Rendering frame 109
[ Info: Rendering frame 110
[ Info: Rendering frame 111
[ Info: Rendering frame 112
[ Info: Rendering frame 113
[ Info: Rendering frame 114
[ Info: Rendering frame 115
[ Info: Rendering frame 116
[ Info: Rendering frame 117
[ Info: Rendering frame 118
[ Info: Rendering frame 119
[ Info: Rendering frame 120
[ Info: Rendering frame 121
[ Info: Rendering frame 122
[ Info: Rendering frame 123
[ Info: Rendering frame 124
[ Info: Rendering frame 125
[ Info: Rendering frame 126
[ Info: Rendering frame 127
[ Info: Rendering frame 128
[ Info: Rendering frame 129
[ Info: Rendering frame 130
[ Info: Rendering frame 131
[ Info: Rendering frame 132
[ Info: Rendering frame 133
[ Info: Rendering frame 134
[ Info: Rendering frame 135
[ Info: Rendering frame 136
[ Info: Rendering frame 137
[ Info: Rendering frame 138
[ Info: Rendering frame 139
[ Info: Rendering frame 140
[ Info: Rendering frame 141
[ Info: Rendering frame 142
[ Info: Rendering frame 143
[ Info: Rendering frame 144
[ Info: Rendering frame 145
[ Info: Rendering frame 146
[ Info: Rendering frame 147
[ Info: Rendering frame 148
[ Info: Rendering frame 149
[ Info: Rendering frame 150
[ Info: Rendering frame 151
[ Info: Rendering frame 152
[ Info: Rendering frame 153
[ Info: Rendering frame 154
[ Info: Rendering frame 155
[ Info: Rendering frame 156
[ Info: Rendering frame 157
[ Info: Rendering frame 158
[ Info: Rendering frame 159
[ Info: Rendering frame 160
[ Info: Rendering frame 161
[ Info: Rendering frame 162
[ Info: Rendering frame 163
[ Info: Rendering frame 164
[ Info: Rendering frame 165
[ Info: Rendering frame 166
[ Info: Rendering frame 167
[ Info: Rendering frame 168
[ Info: Rendering frame 169
[ Info: Rendering frame 170
[ Info: Rendering frame 171
[ Info: Rendering frame 172
[ Info: Rendering frame 173
[ Info: Rendering frame 174
[ Info: Rendering frame 175
[ Info: Rendering frame 176
[ Info: Rendering frame 177
[ Info: Rendering frame 178
[ Info: Rendering frame 179
[ Info: Rendering frame 180
[ Info: Rendering frame 181
[ Info: Rendering frame 182
[ Info: Rendering frame 183
[ Info: Rendering frame 184
[ Info: Rendering frame 185
[ Info: Rendering frame 186
[ Info: Rendering frame 187
[ Info: Rendering frame 188
[ Info: Rendering frame 189
[ Info: Rendering frame 190
[ Info: Rendering frame 191
[ Info: Rendering frame 192
[ Info: Rendering frame 193
[ Info: Rendering frame 194
[ Info: Rendering frame 195
[ Info: Rendering frame 196
[ Info: Rendering frame 197
[ Info: Rendering frame 198
[ Info: Rendering frame 199
[ Info: Rendering frame 200
[ Info: Rendering frame 201
[ Info: Rendering frame 202
[ Info: Rendering frame 203
[ Info: Rendering frame 204
[ Info: Rendering frame 205
[ Info: Rendering frame 206
[ Info: Rendering frame 207
[ Info: Rendering frame 208
[ Info: Rendering frame 209
[ Info: Rendering frame 210
[ Info: Rendering frame 211
[ Info: Rendering frame 212
[ Info: Rendering frame 213
[ Info: Rendering frame 214
[ Info: Rendering frame 215
[ Info: Rendering frame 216
[ Info: Rendering frame 217
[ Info: Rendering frame 218
[ Info: Rendering frame 219
[ Info: Rendering frame 220
[ Info: Rendering frame 221
[ Info: Rendering frame 222
[ Info: Rendering frame 223
[ Info: Rendering frame 224
[ Info: Rendering frame 225
[ Info: Rendering frame 226
[ Info: Rendering frame 227
[ Info: Rendering frame 228
[ Info: Rendering frame 229
[ Info: Rendering frame 230
[ Info: Rendering frame 231
[ Info: Rendering frame 232
[ Info: Rendering frame 233
[ Info: Rendering frame 234
[ Info: Rendering frame 235
[ Info: Rendering frame 236
[ Info: Rendering frame 237
[ Info: Rendering frame 238
[ Info: Rendering frame 239
[ Info: Rendering frame 240
[ Info: Rendering frame 241
[ Info: Rendering frame 242
[ Info: Rendering frame 243
[ Info: Rendering frame 244
[ Info: Rendering frame 245
[ Info: Rendering frame 246
[ Info: Rendering frame 247
[ Info: Rendering frame 248
[ Info: Rendering frame 249
[ Info: Rendering frame 250
[ Info: Rendering frame 251
[ Info: Rendering frame 252
[ Info: Rendering frame 253
[ Info: Rendering frame 254
[ Info: Rendering frame 255
[ Info: Rendering frame 256
[ Info: Rendering frame 257
[ Info: Rendering frame 258
[ Info: Rendering frame 259
[ Info: Rendering frame 260
[ Info: Rendering frame 261
[ Info: Rendering frame 262
[ Info: Rendering frame 263
[ Info: Rendering frame 264
[ Info: Rendering frame 265
[ Info: Rendering frame 266
[ Info: Rendering frame 267
[ Info: Rendering frame 268
[ Info: Rendering frame 269
[ Info: Rendering frame 270
[ Info: Rendering frame 271
[ Info: Rendering frame 272
[ Info: Rendering frame 273
[ Info: Rendering frame 274
[ Info: Rendering frame 275
[ Info: Rendering frame 276
[ Info: Rendering frame 277
[ Info: Rendering frame 278
[ Info: Rendering frame 279
[ Info: Rendering frame 280
[ Info: Rendering frame 281
[ Info: Rendering frame 282
[ Info: Rendering frame 283
[ Info: Rendering frame 284
[ Info: Rendering frame 285
[ Info: Rendering frame 286
[ Info: Rendering frame 287
[ Info: Rendering frame 288
[ Info: Rendering frame 289

The animation shows how the ice is advected and deformed by the moving anticyclonic atmospheric eddy. Linear kinematic features (leads and ridges) form as the ice responds to the complex stress patterns.

This page was generated using Literate.jl.