diff --git a/examples/melting_baroclinicity.jl b/examples/melting_baroclinicity.jl new file mode 100644 index 00000000..fda05588 --- /dev/null +++ b/examples/melting_baroclinicity.jl @@ -0,0 +1,244 @@ +using Oceananigans +using Oceananigans.Architectures: arch_array +using Oceananigans.Fields: ZeroField, ConstantField +using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity +using Oceananigans.Units +using Oceananigans.Utils: prettysummary + +using SeawaterPolynomials: TEOS10EquationOfState, thermal_expansion, haline_contraction + +using ClimaSeaIce +using ClimaSeaIce: melting_temperature +using ClimaSeaIce.ThermalBoundaryConditions: RadiativeEmission, IceWaterThermalEquilibrium + +using Printf +using GLMakie +using Statistics + +include("ice_ocean_model.jl") + +arch = GPU() +Nx = Ny = 256 +Nz = 10 +Lz = 400 +x = y = (-50kilometers, 50kilometers) +halo = (4, 4, 4) +topology = (Periodic, Bounded, Bounded) + +ice_grid = RectilinearGrid(arch; x, y, + size = (Nx, Ny), + topology = (topology[1], topology[2], Flat), + halo = halo[1:2]) + +ocean_grid = RectilinearGrid(arch; topology, halo, x, y, + size = (Nx, Ny, Nz), + z = (-Lz, 0)) + +# Top boundary conditions: +# - outgoing radiative fluxes emitted from surface +# - incoming shortwave radiation starting after 40 days + +ice_ocean_heat_flux = Field{Center, Center, Nothing}(ice_grid) +top_ocean_heat_flux = Qᵀ = Field{Center, Center, Nothing}(ice_grid) +top_salt_flux = Qˢ = Field{Center, Center, Nothing}(ice_grid) +# top_salt_flux = Qˢ = arch_array(arch, zeros(Nx, Ny)) + +# Generate a zero-dimensional grid for a single column slab model + +boundary_conditions = (T = FieldBoundaryConditions(top=FluxBoundaryCondition(Qᵀ)), + S = FieldBoundaryConditions(top=FluxBoundaryCondition(Qˢ))) + +equation_of_state = TEOS10EquationOfState() +buoyancy = SeawaterBuoyancy(; equation_of_state) +horizontal_biharmonic_diffusivity = HorizontalScalarBiharmonicDiffusivity(κ=5e6) + +ocean_model = HydrostaticFreeSurfaceModel(; buoyancy, boundary_conditions, + grid = ocean_grid, + momentum_advection = WENO(), + tracer_advection = WENO(), + #closure = (horizontal_biharmonic_diffusivity, CATKEVerticalDiffusivity()), + closure = CATKEVerticalDiffusivity(), + coriolis = FPlane(f=1.4e-4), + tracers = (:T, :S, :e)) + +Nz = size(ocean_grid, 3) +So = ocean_model.tracers.S +ocean_surface_salinity = view(So, :, :, Nz) +bottom_bc = IceWaterThermalEquilibrium(ConstantField(30)) #ocean_surface_salinity) + +u, v, w = ocean_model.velocities +ocean_surface_velocities = (u = view(u, :, :, Nz), #interior(u, :, :, Nz), + v = view(v, :, :, Nz), #interior(v, :, :, Nz), + w = ZeroField()) + +ice_model = SlabSeaIceModel(ice_grid; + velocities = ocean_surface_velocities, + advection = nothing, #WENO(), + ice_consolidation_thickness = 0.05, + ice_salinity = 4, + internal_thermal_flux = ConductiveFlux(conductivity=2), + #top_thermal_flux = ConstantField(-100), # W m⁻² + top_thermal_flux = ConstantField(0), # W m⁻² + top_thermal_boundary_condition = PrescribedTemperature(0), + bottom_thermal_boundary_condition = bottom_bc, + bottom_thermal_flux = ice_ocean_heat_flux) + +ocean_simulation = Simulation(ocean_model; Δt=20minutes, verbose=false) +ice_simulation = Simulation(ice_model, Δt=20minutes, verbose=false) + +# Initial condition +S₀ = 30 +T₀ = melting_temperature(ice_model.phase_transitions.liquidus, S₀) + 2.0 + +N²S = 1e-6 +β = haline_contraction(T₀, S₀, 0, equation_of_state) +g = ocean_model.buoyancy.model.gravitational_acceleration +dSdz = - g * β * N²S + +uᵢ(x, y, z) = 0.0 +Tᵢ(x, y, z) = T₀ # + 0.1 * randn() +Sᵢ(x, y, z) = S₀ + dSdz * z #+ 0.1 * randn() + +function hᵢ(x, y) + if sqrt(x^2 + y^2) < 20kilometers + #return 1 + 0.1 * rand() + return 2 + else + return 0 + end +end + +set!(ocean_model, u=uᵢ, S=Sᵢ, T=T₀) +set!(ice_model, h=hᵢ) + +coupled_model = IceOceanModel(ice_simulation, ocean_simulation) +coupled_simulation = Simulation(coupled_model, Δt=1minutes, stop_time=20days) + +S = ocean_model.tracers.S +by = - g * β * ∂y(S) + +function progress(sim) + h = sim.model.ice.model.ice_thickness + S = sim.model.ocean.model.tracers.S + T = sim.model.ocean.model.tracers.T + u = sim.model.ocean.model.velocities.u + msg1 = @sprintf("Iter: % 6d, time: % 12s", iteration(sim), prettytime(sim)) + msg2 = @sprintf(", max(h): %.2f", maximum(h)) + msg3 = @sprintf(", min(S): %.2f", minimum(S)) + msg4 = @sprintf(", extrema(T): (%.2f, %.2f)", minimum(T), maximum(T)) + msg5 = @sprintf(", max|∂y b|: %.2e", maximum(abs, by)) + msg6 = @sprintf(", max|u|: %.2e", maximum(abs, u)) + @info msg1 * msg2 * msg3 * msg4 * msg5 * msg6 + return nothing +end + +coupled_simulation.callbacks[:progress] = Callback(progress, IterationInterval(10)) + +h = ice_model.ice_thickness +T = ocean_model.tracers.T +S = ocean_model.tracers.S +u, v, w = ocean_model.velocities +η = ocean_model.free_surface.η + +ht = [] +Tt = [] +Ft = [] +Qt = [] +St = [] +ut = [] +vt = [] +ηt = [] +ζt = [] +tt = [] + +ζ = Field(∂x(v) - ∂y(u)) + +function saveoutput(sim) + compute!(ζ) + hn = Array(interior(h, :, :, 1)) + Fn = Array(interior(Qˢ, :, :, 1)) + Qn = Array(interior(Qᵀ, :, :, 1)) + Tn = Array(interior(T, :, :, Nz)) + Sn = Array(interior(S, :, :, Nz)) + un = Array(interior(u, :, :, Nz)) + vn = Array(interior(v, :, :, Nz)) + ηn = Array(interior(η, :, :, 1)) + ζn = Array(interior(ζ, :, :, Nz)) + push!(ht, hn) + push!(Ft, Fn) + push!(Qt, Qn) + push!(Tt, Tn) + push!(St, Sn) + push!(ut, un) + push!(vt, vn) + push!(ηt, ηn) + push!(ζt, ζn) + push!(tt, time(sim)) +end + +coupled_simulation.callbacks[:output] = Callback(saveoutput, IterationInterval(10)) + +run!(coupled_simulation) + +##### +##### Viz +##### + +set_theme!(Theme(fontsize=24)) + +x = xnodes(ocean_grid, Center()) +y = ynodes(ocean_grid, Center()) + +fig = Figure(resolution=(2400, 700)) + +axh = Axis(fig[1, 1], xlabel="x (km)", ylabel="y (km)", title="Ice thickness") +axT = Axis(fig[1, 2], xlabel="x (km)", ylabel="y (km)", title="Ocean surface temperature") +axS = Axis(fig[1, 3], xlabel="x (km)", ylabel="y (km)", title="Ocean surface salinity") +axZ = Axis(fig[1, 4], xlabel="x (km)", ylabel="y (km)", title="Ocean vorticity") + +Nt = length(tt) +slider = Slider(fig[2, 1:4], range=1:Nt, startvalue=Nt) +n = slider.value + +title = @lift string("Melt-driven baroclinic instability after ", prettytime(tt[$n])) +Label(fig[0, 1:3], title) + +hn = @lift ht[$n] +Fn = @lift Ft[$n] +Tn = @lift Tt[$n] +Sn = @lift St[$n] +un = @lift ut[$n] +vn = @lift vt[$n] +ηn = @lift ηt[$n] +ζn = @lift ζt[$n] +Un = @lift mean(ut[$n], dims=1)[:] + +x = x ./ 1e3 +y = y ./ 1e3 + +Stop = view(S, :, :, Nz) +Smax = maximum(Stop) +Smin = minimum(Stop) + +compute!(ζ) +ζtop = view(ζ, :, :, Nz) +ζmax = maximum(abs, ζtop) +ζlim = 2e-4 #ζmax / 2 + +heatmap!(axh, x, y, hn, colorrange=(0, 1), colormap=:grays) +heatmap!(axT, x, y, Tn, colormap=:thermal) +heatmap!(axS, x, y, Sn, colorrange = (29, 30), colormap=:haline) +heatmap!(axZ, x, y, ζn, colorrange=(-ζlim, ζlim), colormap=:redblue) + +#heatmap!(axZ, x, y, Tn, colormap=:thermal) +#heatmap!(axF, x, y, Fn) + +display(fig) + +#= +record(fig, "salty_baroclinic_ice_cube.mp4", 1:Nt, framerate=48) do nn + @info string(nn) + n[] = nn +end +=# + diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index 64a8aaf9..903436ab 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -105,5 +105,6 @@ include("DataWrangling.jl") include("Diagnostics.jl") include("NearGlobalSimulations/NearGlobalSimulations.jl") include("IdealizedSimulations/IdealizedSimulations.jl") +include("IceOceanModel/IceOceanModel.jl") end # module diff --git a/src/IceOceanModel/IceOceanModel.jl b/src/IceOceanModel/IceOceanModel.jl new file mode 100644 index 00000000..d4b9696e --- /dev/null +++ b/src/IceOceanModel/IceOceanModel.jl @@ -0,0 +1,348 @@ +module IceOceanModel + +using Oceananigans.Operators + +using Oceananigans.Architectures: architecture +using Oceananigans.BoundaryConditions: fill_halo_regions! +using Oceananigans.Models: AbstractModel +using Oceananigans.TimeSteppers: tick! +using Oceananigans.Utils: launch! + +using KernelAbstractions: @kernel, @index +using KernelAbstractions.Extras.LoopInfo: @unroll + +# Simulations interface +import Oceananigans: fields, prognostic_fields +import Oceananigans.Fields: set! +import Oceananigans.Models: timestepper, NaNChecker, default_nan_checker +import Oceananigans.OutputWriters: default_included_properties +import Oceananigans.Simulations: reset!, initialize!, iteration +import Oceananigans.TimeSteppers: time_step!, update_state!, time +import Oceananigans.Utils: prettytime + +struct IceOceanModel{FT, C, G, I, O, S, PI, PC} <: AbstractModel{Nothing} + clock :: C + grid :: G # TODO: make it so simulation does not require this + ice :: I + previous_ice_thickness :: PI + previous_ice_concentration :: PC + ocean :: O + solar_insolation :: S + ocean_density :: FT + ocean_heat_capacity :: FT + ocean_emissivity :: FT + stefan_boltzmann_constant :: FT + reference_temperature :: FT +end + +const IOM = IceOceanModel + +Base.summary(::IOM) = "IceOceanModel" +prettytime(model::IOM) = prettytime(model.clock.time) +iteration(model::IOM) = model.clock.iteration +timestepper(::IOM) = nothing +reset!(::IOM) = nothing +initialize!(::IOM) = nothing +default_included_properties(::IOM) = tuple() +update_state!(::IOM) = nothing +prognostic_fields(cm::IOM) = nothing +fields(::IOM) = NamedTuple() + +function IceOceanModel(ice, ocean; clock = Clock{Float64}(0, 0, 1)) + + previous_ice_thickness = deepcopy(ice.model.ice_thickness) + previous_ice_concentration = deepcopy(ice.model.ice_concentration) + + grid = ocean.model.grid + ice_ocean_thermal_flux = Field{Center, Center, Nothing}(grid) + ice_ocean_salt_flux = Field{Center, Center, Nothing}(grid) + solar_insolation = Field{Center, Center, Nothing}(grid) + + ocean_density = 1024 + ocean_heat_capacity = 3991 + ocean_emissivity = 1 + reference_temperature = 273.15 + stefan_boltzmann_constant = 5.67e-8 + + # How would we ensure consistency? + try + if ice.model.external_thermal_fluxes.top isa RadiativeEmission + radiation = ice.model.external_thermal_fluxes.top + else + radiation = filter(flux isa RadiativeEmission, ice.model.external_thermal_fluxes.top) |> first + end + + stefan_boltzmann_constant = radiation.stefan_boltzmann_constant + reference_temperature = radiation.reference_temperature + catch + end + + FT = eltype(ocean.model.grid) + + return IceOceanModel(clock, + ocean.model.grid, + ice, + previous_ice_thickness, + previous_ice_concentration, + ocean, + solar_insolation, + convert(FT, ocean_density), + convert(FT, ocean_heat_capacity), + convert(FT, ocean_emissivity), + convert(FT, stefan_boltzmann_constant), + convert(FT, reference_temperature)) +end + +time(coupled_model::IceOceanModel) = coupled_model.clock.time + +function compute_air_sea_flux!(coupled_model) + ocean = coupled_model.ocean + ice = coupled_model.ice + + T = ocean.model.tracers.T + Nx, Ny, Nz = size(ocean.model.grid) + + grid = ocean.model.grid + arch = architecture(grid) + + σ = coupled_model.stefan_boltzmann_constant + ρₒ = coupled_model.ocean_density + cₒ = coupled_model.ocean_heat_capacity + Tᵣ = coupled_model.reference_temperature + Qᵀ = T.boundary_conditions.top.condition + Tₒ = ocean.model.tracers.T + hᵢ = ice.model.ice_thickness + ℵᵢ = ice.model.ice_concentration + I₀ = coupled_model.solar_insolation + + launch!(arch, grid, :xy, _compute_air_sea_flux!, + Qᵀ, grid, Tₒ, hᵢ, ℵᵢ, I₀, + σ, ρₒ, cₒ, Tᵣ) +end + +@kernel function _compute_air_sea_flux!(temperature_flux, + grid, + ocean_temperature, + ice_thickness, + ice_concentration, + solar_insolation, + σ, ρₒ, cₒ, Tᵣ) + + i, j = @index(Global, NTuple) + + Nz = size(grid, 3) + + @inbounds begin + T₀ = ocean_temperature[i, j, Nz] # at the surface + h = ice_thickness[i, j, 1] + ℵ = ice_concentration[i, j, 1] + I₀ = solar_insolation[i, j, 1] + end + + # Radiation model + ϵ = 1 # ocean emissivity + ΣQᵀ = ϵ * σ * (T₀ + Tᵣ)^4 / (ρₒ * cₒ) + + # Also add solar insolation + ΣQᵀ += I₀ / (ρₒ * cₒ) + + # Set the surface flux only if ice-free + Qᵀ = temperature_flux + + # @inbounds Qᵀ[i, j, 1] = 0 #(1 - ℵ) * ΣQᵀ +end + +function time_step!(coupled_model::IceOceanModel, Δt; callbacks=nothing) + ocean = coupled_model.ocean + ice = coupled_model.ice + liquidus = ice.model.phase_transitions.liquidus + grid = ocean.model.grid + ice.Δt = Δt + ocean.Δt = Δt + + fill_halo_regions!(h) + + # Initialization + if coupled_model.clock.iteration == 0 + h⁻ = coupled_model.previous_ice_thickness + hⁿ = coupled_model.ice.model.ice_thickness + parent(h⁻) .= parent(hⁿ) + end + + time_step!(ice) + + # TODO: put this in update_state! + compute_ice_ocean_salinity_flux!(coupled_model) + ice_ocean_latent_heat!(coupled_model) + #compute_solar_insolation!(coupled_model) + #compute_air_sea_flux!(coupled_model) + + time_step!(ocean) + + # TODO: + # - Store fractional ice-free / ice-covered _time_ for more + # accurate flux computation? + # - Or, input "excess heat flux" into ocean after the ice melts + # - Currently, non-conservative for heat due bc we don't account for excess + + # TODO after ice time-step: + # - Adjust ocean temperature if the ice completely melts? + + tick!(coupled_model.clock, Δt) + + return nothing +end + +function compute_ice_ocean_salinity_flux!(coupled_model) + # Compute salinity increment due to changes in ice thickness + + ice = coupled_model.ice + ocean = coupled_model.ocean + grid = ocean.model.grid + arch = architecture(grid) + Qˢ = ocean.model.tracers.S.boundary_conditions.top.condition + Sₒ = ocean.model.tracers.S + Sᵢ = ice.model.ice_salinity + Δt = ocean.Δt + hⁿ = ice.model.ice_thickness + h⁻ = coupled_model.previous_ice_thickness + + launch!(arch, grid, :xy, _compute_ice_ocean_salinity_flux!, + Qˢ, grid, hⁿ, h⁻, Sᵢ, Sₒ, Δt) + + return nothing +end + + +@kernel function _compute_ice_ocean_salinity_flux!(ice_ocean_salinity_flux, + grid, + ice_thickness, + previous_ice_thickness, + ice_salinity, + ocean_salinity, + Δt) + i, j = @index(Global, NTuple) + + Nz = size(grid, 3) + + hⁿ = ice_thickness + h⁻ = previous_ice_thickness + Qˢ = ice_ocean_salinity_flux + Sᵢ = ice_salinity + Sₒ = ocean_salinity + + @inbounds begin + # Thickness of surface grid cell + Δh = hⁿ[i, j, 1] - h⁻[i, j, 1] + + # Update surface salinity flux. + # Note: the Δt below is the ocean time-step, eg. + # ΔS = ⋯ - ∮ Qˢ dt ≈ ⋯ - Δtₒ * Qˢ + Qˢ[i, j, 1] = Δh / Δt * (Sᵢ[i, j, 1] - Sₒ[i, j, Nz]) + + # Update previous ice thickness + h⁻[i, j, 1] = hⁿ[i, j, 1] + end +end + +function ice_ocean_latent_heat!(coupled_model) + ocean = coupled_model.ocean + ice = coupled_model.ice + ρₒ = coupled_model.ocean_density + cₒ = coupled_model.ocean_heat_capacity + Qₒ = ice.model.external_thermal_fluxes.bottom + Tₒ = ocean.model.tracers.T + Sₒ = ocean.model.tracers.S + Δt = ocean.Δt + hᵢ = ice.model.ice_thickness + + liquidus = ice.model.phase_transitions.liquidus + grid = ocean.model.grid + arch = architecture(grid) + + launch!(arch, grid, :xy, _compute_ice_ocean_latent_heat!, + Qₒ, grid, hᵢ, Tₒ, Sₒ, liquidus, ρₒ, cₒ, Δt) + + return nothing +end + +@kernel function _compute_ice_ocean_latent_heat!(latent_heat, + grid, + ice_thickness, + ocean_temperature, + ocean_salinity, + liquidus, + ρₒ, cₒ, Δt) + + i, j = @index(Global, NTuple) + + Nz = size(grid, 3) + Qₒ = latent_heat + hᵢ = ice_thickness + Tₒ = ocean_temperature + Sₒ = ocean_salinity + + δQ = zero(grid) + icy_cell = @inbounds hᵢ[i, j, 1] > 0 # make ice bath approximation then + + @unroll for k = Nz:-1:1 + @inbounds begin + # Various quantities + Δz = Δzᶜᶜᶜ(i, j, k, grid) + Tᴺ = Tₒ[i, j, k] + Sᴺ = Sₒ[i, j, k] + end + + # Melting / freezing temperature at the surface of the ocean + Tₘ = melting_temperature(liquidus, Sᴺ) + + # Conditions for non-zero ice-ocean flux: + # - the ocean is below the freezing temperature, causing formation of ice. + freezing = Tᴺ < Tₘ + + # - We are at the surface and the cell is covered by ice. + icy_surface_cell = (k == Nz) & icy_cell + + # When there is a non-zero ice-ocean flux, we will instantaneously adjust the + # temperature of the grid cells accordingly. + adjust_temperature = freezing | icy_surface_cell + + # Compute change in ocean thermal energy. + # + # - When Tᴺ < Tₘ, we heat the ocean back to melting temperature by extracting heat from the ice, + # assuming that the heat flux (which is carried by nascent ice crystals called frazil ice) floats + # instantaneously to the surface. + # + # - When Tᴺ > Tₘ and we are in a surface cell covered by ice, we assume equilibrium + # and cool the ocean by injecting excess heat into the ice. + # + δEₒ = adjust_temperature * ρₒ * cₒ * (Tₘ - Tᴺ) + + # Perform temperature adjustment + @inline Tₒ[i, j, k] = ifelse(adjust_temperature, Tₘ, Tᴺ) + + # Compute the heat flux from ocean into ice. + # + # A positive value δQ > 0 implies that the ocean is cooled; ie heat + # is fluxing upwards, into the ice. This occurs when applying the + # ice bath equilibrium condition to cool down a warm ocean (δEₒ < 0). + # + # A negative value δQ < 0 implies that heat is fluxed from the ice into + # the ocean, cooling the ice and heating the ocean (δEₒ > 0). This occurs when + # frazil ice is formed within the ocean. + + δQ -= δEₒ * Δz / Δt + end + + # Store ice-ocean flux + @inbounds Qₒ[i, j, 1] = δQ +end + +# Check for NaNs in the first prognostic field (generalizes to prescribed velocitries). +function default_nan_checker(model::IceOceanModel) + u_ocean = model.ocean.model.velocities.u + nan_checker = NaNChecker((; u_ocean)) + return nan_checker +end + +end # module