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 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
└── 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 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=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
                    timestepper = :SplitRungeKutta3,
                    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

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
├── 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

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",
                                                 including = [:grid],
                                                 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(100))

run!(simulation)

[ Info: Initializing simulation...
[ Info: Time: 0 seconds, Iteration 0, Δt 2 minutes, max(vel): (0.00e+00, 0.00e+00), max(h): 0.31, min(ℵ): 1.00, wtime: 7.225 seconds
[ Info:     ... simulation initialization complete (32.130 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (7.651 seconds).
[ Info: Time: 3.333 hours, Iteration 100, Δt 2 minutes, max(vel): (1.55e-01, 1.55e-01), max(h): 0.31, min(ℵ): 0.98, wtime: 2.176 minutes
[ Info: Time: 6.667 hours, Iteration 200, Δt 2 minutes, max(vel): (1.55e-01, 1.55e-01), max(h): 0.31, min(ℵ): 0.97, wtime: 1.620 minutes
[ Info: Time: 10 hours, Iteration 300, Δt 2 minutes, max(vel): (1.55e-01, 1.55e-01), max(h): 0.32, min(ℵ): 0.95, wtime: 1.626 minutes
[ Info: Time: 13.333 hours, Iteration 400, Δt 2 minutes, max(vel): (1.55e-01, 1.55e-01), max(h): 0.32, min(ℵ): 0.92, wtime: 1.631 minutes
[ Info: Time: 16.667 hours, Iteration 500, Δt 2 minutes, max(vel): (1.55e-01, 1.55e-01), max(h): 0.32, min(ℵ): 0.90, wtime: 1.628 minutes
[ Info: Time: 20 hours, Iteration 600, Δt 2 minutes, max(vel): (1.55e-01, 1.56e-01), max(h): 0.32, min(ℵ): 0.85, wtime: 1.630 minutes
[ Info: Time: 23.333 hours, Iteration 700, Δt 2 minutes, max(vel): (1.55e-01, 1.56e-01), max(h): 0.32, min(ℵ): 0.78, wtime: 1.622 minutes
[ Info: Time: 1.111 days, Iteration 800, Δt 2 minutes, max(vel): (1.55e-01, 1.57e-01), max(h): 0.33, min(ℵ): 0.70, wtime: 1.627 minutes
[ Info: Time: 1.250 days, Iteration 900, Δt 2 minutes, max(vel): (1.55e-01, 1.59e-01), max(h): 0.33, min(ℵ): 0.61, wtime: 1.622 minutes
[ Info: Time: 1.389 days, Iteration 1000, Δt 2 minutes, max(vel): (1.55e-01, 1.61e-01), max(h): 0.34, min(ℵ): 0.51, wtime: 1.625 minutes
[ Info: Time: 1.528 days, Iteration 1100, Δt 2 minutes, max(vel): (1.55e-01, 1.66e-01), max(h): 0.35, min(ℵ): 0.40, wtime: 1.630 minutes
[ Info: Time: 1.667 days, Iteration 1200, Δt 2 minutes, max(vel): (1.55e-01, 1.71e-01), max(h): 0.36, min(ℵ): 0.29, wtime: 1.632 minutes
[ Info: Time: 1.806 days, Iteration 1300, Δt 2 minutes, max(vel): (1.55e-01, 1.75e-01), max(h): 0.37, min(ℵ): 0.17, wtime: 1.630 minutes
[ Info: Time: 1.944 days, Iteration 1400, Δt 2 minutes, max(vel): (1.56e-01, 1.79e-01), max(h): 0.38, min(ℵ): 0.06, wtime: 1.628 minutes
[ Info: Simulation is stopping after running for 24.045 minutes.
[ Info: Simulation time 2 days equals or exceeds stop time 2 days.

Load output

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, Bounded, Bounded, Flat} on 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()
ax1 = Axis(fig[1, 1], title = "sea ice thickness")
ax2 = Axis(fig[1, 2], title = "sea ice concentration")
ax3 = Axis(fig[2, 1], title = "zonal velocity")
ax4 = Axis(fig[2, 2], title = "meridional velocity")

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))
heatmap!(ax4, vi, colorrange = (-0.1, 0.1))

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

This page was generated using Literate.jl.