NumericalEarth · simone-silvestri · Apr 15, 2026 · Apr 15, 2026 · Apr 15, 2026 · Apr 20, 2026
diff --git a/src/Bathymetry/regrid_bathymetry.jl b/src/Bathymetry/regrid_bathymetry.jl
@@ -341,7 +341,8 @@ end
 
     # Fix active cells to be at least `-minimum_depth`.
     active = z < 0 # it's a wet cell
-    z = ifelse(active, min(z, -minimum_depth), z)
+    above_minimum_depth = z > -minimum_depth
+    z = ifelse(active, ifelse(above_minimum_depth, zero(z), z), z)
 
     @inbounds target_z[i, j, 1] = z
 end

diff --git a/src/DataWrangling/restoring.jl b/src/DataWrangling/restoring.jl
@@ -14,6 +14,7 @@ using Dates: Second
 
 import NumericalEarth: stateindex
 import Oceananigans.Forcings: materialize_forcing
+import Oceananigans.OutputReaders: extract_field_time_series
 
 # Variable names for restorable data
 struct Temperature end
@@ -236,6 +237,7 @@ function Base.show(io::IO, dsr::DatasetRestoring)
 end
 
 materialize_forcing(forcing::DatasetRestoring, field, field_name, model_field_names) = forcing
+extract_field_time_series(forcing::DatasetRestoring) = forcing.field_time_series
 
 """
     SurfaceFluxRestoring(dataset_restoring::DatasetRestoring)

diff --git a/src/Diagnostics/mixed_layer_depth.jl b/src/Diagnostics/mixed_layer_depth.jl
@@ -32,7 +32,6 @@ end
 
 function compute!(mld::MixedLayerDepthField, time=nothing)
     compute_mixed_layer_depth!(mld)
-    #@apply_regionally compute_mixed_layer_depth!(mld)
     fill_halo_regions!(mld)
     return mld
 end

diff --git a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_fluxes.jl
@@ -175,12 +175,12 @@ end
     qᶠ = δ𝒬ᶠʳᶻ / ℰ
 
     @inbounds begin
-        Tᴺ = Tᵒᶜ[i, j, Nz]               
-        Sᴺ = Sᵒᶜ[i, j, Nz]               
+        Tᴺ  = Tᵒᶜ[i, j, Nz]               
+        Sᴺ  = Sᵒᶜ[i, j, Nz]               
         Sˢⁱ = ice_salinity[i, j, 1]      
         hˢⁱ = ice_thickness[i, j, 1]     
-        ℵᵢ = ice_concentration[i, j, 1] 
-        hc = ice_consolidation_thickness[i, j, 1] 
+        ℵᵢ  = ice_concentration[i, j, 1] 
+        hc  = ice_consolidation_thickness[i, j, 1] 
     end
 
     # Extract internal temperature (for ConductiveFluxTEF, zero otherwise)
@@ -198,8 +198,8 @@ end
     # =============================================
     # Returns interfacial heat flux, melt rate qᵐ, and interface T, S
     𝒬ⁱᵒ, qᵐ, Tᵦ, Sᵦ = compute_interface_heat_flux(flux_formulation,
-                                                     ocean_surface_state, ice_state,
-                                                     liquidus, ocean_properties, ℰ, u★)
+                                                   ocean_surface_state, ice_state,
+                                                   liquidus, ocean_properties, ℰ, u★)
 
     # Store interface values and heat flux
     @inbounds 𝒬ⁱⁿᵗ[i, j, 1] = 𝒬ⁱᵒ

diff --git a/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_heat_flux_formulations.jl b/src/EarthSystemModels/InterfaceComputations/sea_ice_ocean_heat_flux_formulations.jl
@@ -113,8 +113,8 @@ References
 
 - [holland1999modeling](@citet): Holland, D. M., & Jenkins, A. (1999). Modeling thermodynamic ice–ocean interactions
   at the base of an ice shelf. *Journal of Physical Oceanography*, 29(8), 1787-1800.
-- [hieronymus2021comparison](@citet): Hieronymus, M., et al. (2021). A comparison of ocean-ice flux parametrizations.
-  *Geosci. Model Dev.*, 14, 4891-4908.
+- [shi2021sensitivity](@citet): Shi, X., Notz, D., Liu, J., Yang, H., & Lohmann, G. (2021). Sensitivity of Northern
+  Hemisphere climate to ice-ocean interface heat flux parameterizations. *Geosci. Model Dev.*, 14, 4891-4908.
 """
 struct ThreeEquationHeatFlux{F, T, FT, U}
     conductive_flux :: F
@@ -139,7 +139,7 @@ Adapt.adapt_structure(to, f::ThreeEquationHeatFlux) =
 
 Construct a `ThreeEquationHeatFlux` with the specified parameters.
 
-Default values follow [hieronymus2021comparison](@citet) with ``R = \\alpha_h / \\alpha_s = 35``.
+Default values follow [shi2021sensitivity](@citet) with ``R = \\alpha_h / \\alpha_s = 35``.
 
 Keyword Arguments
 =================

diff --git a/src/EarthSystemModels/earth_system_model.jl b/src/EarthSystemModels/earth_system_model.jl
@@ -44,13 +44,12 @@ function Base.show(io::IO, cm::ESM)
     return nothing
 end
 
-# Assumption: We have an ocean!
 architecture(model::ESM)           = model.architecture
 Base.eltype(model::ESM)            = Base.eltype(model.interfaces.exchanger.grid)
 prettytime(model::ESM)             = prettytime(model.clock.time)
 iteration(model::ESM)              = model.clock.iteration
 timestepper(::ESM)                 = nothing
-default_included_properties(::ESM) = tuple()
+default_included_properties(::ESM) = []
 prognostic_fields(cm::ESM)         = nothing
 fields(::ESM)                      = NamedTuple()
 default_clock(TT)                   = Oceananigans.TimeSteppers.Clock{TT}(0, 0, 1)

diff --git a/src/EarthSystemModels/time_step_earth_system_model.jl b/src/EarthSystemModels/time_step_earth_system_model.jl
@@ -15,13 +15,8 @@ function time_step!(coupled_model::EarthSystemModel, Δt; callbacks=[])
     atmosphere = coupled_model.atmosphere
 
     # Eventually, split out into OceanOnlyModel
-    !isnothing(sea_ice) && time_step!(sea_ice, Δt)
-
-    # TODO after ice time-step:
-    #  - Adjust ocean heat flux if the ice completely melts?
-    !isnothing(ocean) && time_step!(ocean, Δt)
-
-    # Time step the atmosphere
+    !isnothing(sea_ice)    && time_step!(sea_ice, Δt)
+    !isnothing(ocean)      && time_step!(ocean, Δt)    
     !isnothing(atmosphere) && time_step!(atmosphere, Δt)
 
     # TODO: