Sea ice advected by an atmospheric anticyclone

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

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

arch = CPU()

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

30.0

2 km domain

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

#####
##### Value boundary conditions for velocities
#####

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

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

#####
##### Ocean sea-ice stress
#####

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)

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

####
#### Atmosphere - sea ice stress
####

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

128×129×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on Oceananigans.Architectures.CPU
├── grid: 128×128×1 RectilinearGrid{Float64, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded, Oceananigans.Grids.Flat} on Oceananigans.Architectures.CPU with 7×7×0 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: ZeroFlux, east: ZeroFlux, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: Nothing
└── 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.0, min=0.0, mean=0.0

Atmospheric velocities corresponding to an anticyclonic eddy moving north-east

@inline center(t) = 256kilometers + 51.2kilometers * t / 86400
@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) * (cosd(72) * (y - center(t)) - sind(72) * (x - center(t))) / 1000

va_time (generic function with 1 method)

Initialize the stress at time t = 0

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)

#####
##### Numerical details
#####

128×129×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on Oceananigans.Architectures.CPU
├── grid: 128×128×1 RectilinearGrid{Float64, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded, Oceananigans.Grids.Flat} on Oceananigans.Architectures.CPU with 7×7×0 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: ZeroFlux, east: ZeroFlux, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: Nothing
├── operand: BinaryOperation 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=3.36248e-8

We use an elasto-visco-plastic rheology and WENO seventh order for advection of h and ℵ

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 ice_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, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded, Oceananigans.Grids.Flat} on Oceananigans.Architectures.CPU with 7×7×0 halo
├── ice_thermodynamics: Nothing
├── advection: WENO{4, Float64, Float32}(order=7)
└── external_heat_fluxes: 
    ├── top: Nothing
    └── bottom: 0

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

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

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

#####
##### Setup the simulation
#####

run the model for 2 days

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

Simulation of SeaIceModel
├── Next time step: 2 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per 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 => 3
│   ├── stop_iteration_exceeded => -
│   └── wall_time_limit_exceeded => }
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entries

Remember to evolve the wind stress field in time!

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

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)

wall_time = [time_ns()]

function progress(sim)
    h = sim.model.ice_thickness
    ℵ = sim.model.ice_concentration
    u = sim.model.velocities.u
    v = sim.model.velocities.v

    hmax = maximum(interior(h))
    ℵmin = minimum(interior(ℵ))
    umax = maximum(interior(u)), maximum(interior(v))
    step_time = 1e-9 * (time_ns() - wall_time[1])

    @info @sprintf("Time: %s, Iteration %d, Δt %s, max(vel): (%.2e, %.2e), max(h): %.2f, min(ℵ): %.2f, wtime: %s \n",
                   prettytime(sim.model.clock.time),
                   sim.model.clock.iteration,
                   prettytime(sim.Δt),
                   umax..., hmax, ℵmin, prettytime(step_time))

     wall_time[1] = time_ns()
end

simulation.callbacks[:progress] = Callback(progress, IterationInterval(5))

run!(simulation)

using CairoMakie

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", "ℵ")

128×128×1×289 FieldTimeSeries{Oceananigans.OutputReaders.InMemory} located at (Center, Center, ⋅) of ℵ at sea_ice_advected_by_anticyclone.jld2
├── grid: 128×128×1 RectilinearGrid{Float64, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded, Oceananigans.Grids.Flat} on Oceananigans.Architectures.CPU with 7×7×0 halo
├── indices: (:, :, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: sea_ice_advected_by_anticyclone.jld2
├── name: ℵ
└── data: 142×142×1×289 OffsetArray(::Array{Float64, 4}, -6:135, -6:135, 1:1, 1:289) with eltype Float64 with indices -6:135×-6:135×1:1×1:289
    └── max=1.0, min=0.0, mean=0.83238

Visualize!

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()
ax = Axis(fig[1, 1], title = "sea ice thickness")
heatmap!(ax, hi, colormap = :magma, colorrange = (0.23, 0.37))

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

ax = Axis(fig[2, 1], title = "zonal velocity")
heatmap!(ax, ui, colorrange = (-0.1, 0.1))

ax = Axis(fig[2, 2], title = "meridional velocity")
heatmap!(ax, vi, colorrange = (-0.1, 0.1))

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

[ Info: doing iter 1
[ Info: doing iter 2
[ Info: doing iter 3
[ Info: doing iter 4
[ Info: doing iter 5
[ Info: doing iter 6
[ Info: doing iter 7
[ Info: doing iter 8
[ Info: doing iter 9
[ Info: doing iter 10
[ Info: doing iter 11
[ Info: doing iter 12
[ Info: doing iter 13
[ Info: doing iter 14
[ Info: doing iter 15
[ Info: doing iter 16
[ Info: doing iter 17
[ Info: doing iter 18
[ Info: doing iter 19
[ Info: doing iter 20
[ Info: doing iter 21
[ Info: doing iter 22
[ Info: doing iter 23
[ Info: doing iter 24
[ Info: doing iter 25
[ Info: doing iter 26
[ Info: doing iter 27
[ Info: doing iter 28
[ Info: doing iter 29
[ Info: doing iter 30
[ Info: doing iter 31
[ Info: doing iter 32
[ Info: doing iter 33
[ Info: doing iter 34
[ Info: doing iter 35
[ Info: doing iter 36
[ Info: doing iter 37
[ Info: doing iter 38
[ Info: doing iter 39
[ Info: doing iter 40
[ Info: doing iter 41
[ Info: doing iter 42
[ Info: doing iter 43
[ Info: doing iter 44
[ Info: doing iter 45
[ Info: doing iter 46
[ Info: doing iter 47
[ Info: doing iter 48
[ Info: doing iter 49
[ Info: doing iter 50
[ Info: doing iter 51
[ Info: doing iter 52
[ Info: doing iter 53
[ Info: doing iter 54
[ Info: doing iter 55
[ Info: doing iter 56
[ Info: doing iter 57
[ Info: doing iter 58
[ Info: doing iter 59
[ Info: doing iter 60
[ Info: doing iter 61
[ Info: doing iter 62
[ Info: doing iter 63
[ Info: doing iter 64
[ Info: doing iter 65
[ Info: doing iter 66
[ Info: doing iter 67
[ Info: doing iter 68
[ Info: doing iter 69
[ Info: doing iter 70
[ Info: doing iter 71
[ Info: doing iter 72
[ Info: doing iter 73
[ Info: doing iter 74
[ Info: doing iter 75
[ Info: doing iter 76
[ Info: doing iter 77
[ Info: doing iter 78
[ Info: doing iter 79
[ Info: doing iter 80
[ Info: doing iter 81
[ Info: doing iter 82
[ Info: doing iter 83
[ Info: doing iter 84
[ Info: doing iter 85
[ Info: doing iter 86
[ Info: doing iter 87
[ Info: doing iter 88
[ Info: doing iter 89
[ Info: doing iter 90
[ Info: doing iter 91
[ Info: doing iter 92
[ Info: doing iter 93
[ Info: doing iter 94
[ Info: doing iter 95
[ Info: doing iter 96
[ Info: doing iter 97
[ Info: doing iter 98
[ Info: doing iter 99
[ Info: doing iter 100
[ Info: doing iter 101
[ Info: doing iter 102
[ Info: doing iter 103
[ Info: doing iter 104
[ Info: doing iter 105
[ Info: doing iter 106
[ Info: doing iter 107
[ Info: doing iter 108
[ Info: doing iter 109
[ Info: doing iter 110
[ Info: doing iter 111
[ Info: doing iter 112
[ Info: doing iter 113
[ Info: doing iter 114
[ Info: doing iter 115
[ Info: doing iter 116
[ Info: doing iter 117
[ Info: doing iter 118
[ Info: doing iter 119
[ Info: doing iter 120
[ Info: doing iter 121
[ Info: doing iter 122
[ Info: doing iter 123
[ Info: doing iter 124
[ Info: doing iter 125
[ Info: doing iter 126
[ Info: doing iter 127
[ Info: doing iter 128
[ Info: doing iter 129
[ Info: doing iter 130
[ Info: doing iter 131
[ Info: doing iter 132
[ Info: doing iter 133
[ Info: doing iter 134
[ Info: doing iter 135
[ Info: doing iter 136
[ Info: doing iter 137
[ Info: doing iter 138
[ Info: doing iter 139
[ Info: doing iter 140
[ Info: doing iter 141
[ Info: doing iter 142
[ Info: doing iter 143
[ Info: doing iter 144
[ Info: doing iter 145
[ Info: doing iter 146
[ Info: doing iter 147
[ Info: doing iter 148
[ Info: doing iter 149
[ Info: doing iter 150
[ Info: doing iter 151
[ Info: doing iter 152
[ Info: doing iter 153
[ Info: doing iter 154
[ Info: doing iter 155
[ Info: doing iter 156
[ Info: doing iter 157
[ Info: doing iter 158
[ Info: doing iter 159
[ Info: doing iter 160
[ Info: doing iter 161
[ Info: doing iter 162
[ Info: doing iter 163
[ Info: doing iter 164
[ Info: doing iter 165
[ Info: doing iter 166
[ Info: doing iter 167
[ Info: doing iter 168
[ Info: doing iter 169
[ Info: doing iter 170
[ Info: doing iter 171
[ Info: doing iter 172
[ Info: doing iter 173
[ Info: doing iter 174
[ Info: doing iter 175
[ Info: doing iter 176
[ Info: doing iter 177
[ Info: doing iter 178
[ Info: doing iter 179
[ Info: doing iter 180
[ Info: doing iter 181
[ Info: doing iter 182
[ Info: doing iter 183
[ Info: doing iter 184
[ Info: doing iter 185
[ Info: doing iter 186
[ Info: doing iter 187
[ Info: doing iter 188
[ Info: doing iter 189
[ Info: doing iter 190
[ Info: doing iter 191
[ Info: doing iter 192
[ Info: doing iter 193
[ Info: doing iter 194
[ Info: doing iter 195
[ Info: doing iter 196
[ Info: doing iter 197
[ Info: doing iter 198
[ Info: doing iter 199
[ Info: doing iter 200
[ Info: doing iter 201
[ Info: doing iter 202
[ Info: doing iter 203
[ Info: doing iter 204
[ Info: doing iter 205
[ Info: doing iter 206
[ Info: doing iter 207
[ Info: doing iter 208
[ Info: doing iter 209
[ Info: doing iter 210
[ Info: doing iter 211
[ Info: doing iter 212
[ Info: doing iter 213
[ Info: doing iter 214
[ Info: doing iter 215
[ Info: doing iter 216
[ Info: doing iter 217
[ Info: doing iter 218
[ Info: doing iter 219
[ Info: doing iter 220
[ Info: doing iter 221
[ Info: doing iter 222
[ Info: doing iter 223
[ Info: doing iter 224
[ Info: doing iter 225
[ Info: doing iter 226
[ Info: doing iter 227
[ Info: doing iter 228
[ Info: doing iter 229
[ Info: doing iter 230
[ Info: doing iter 231
[ Info: doing iter 232
[ Info: doing iter 233
[ Info: doing iter 234
[ Info: doing iter 235
[ Info: doing iter 236
[ Info: doing iter 237
[ Info: doing iter 238
[ Info: doing iter 239
[ Info: doing iter 240
[ Info: doing iter 241
[ Info: doing iter 242
[ Info: doing iter 243
[ Info: doing iter 244
[ Info: doing iter 245
[ Info: doing iter 246
[ Info: doing iter 247
[ Info: doing iter 248
[ Info: doing iter 249
[ Info: doing iter 250
[ Info: doing iter 251
[ Info: doing iter 252
[ Info: doing iter 253
[ Info: doing iter 254
[ Info: doing iter 255
[ Info: doing iter 256
[ Info: doing iter 257
[ Info: doing iter 258
[ Info: doing iter 259
[ Info: doing iter 260
[ Info: doing iter 261
[ Info: doing iter 262
[ Info: doing iter 263
[ Info: doing iter 264
[ Info: doing iter 265
[ Info: doing iter 266
[ Info: doing iter 267
[ Info: doing iter 268
[ Info: doing iter 269
[ Info: doing iter 270
[ Info: doing iter 271
[ Info: doing iter 272
[ Info: doing iter 273
[ Info: doing iter 274
[ Info: doing iter 275
[ Info: doing iter 276
[ Info: doing iter 277
[ Info: doing iter 278
[ Info: doing iter 279
[ Info: doing iter 280
[ Info: doing iter 281
[ Info: doing iter 282
[ Info: doing iter 283
[ Info: doing iter 284
[ Info: doing iter 285
[ Info: doing iter 286
[ Info: doing iter 287
[ Info: doing iter 288
[ Info: doing iter 289

This page was generated using Literate.jl.