Implement a generalized column thermodynamics model

[glwagner]

Summary

Adds column-energy sea-ice thermodynamics with the fixed-salinity Bitz-Lipscomb BL99 reproduction path, evolving-salinity mushy-column configuration, and CICE/Icepack single-column validation evidence.

This PR now includes:

• ColumnEnergyThermodynamics field containers, presets, checkpoint helpers, Oceananigans integration, and a focused column_energy test group.
• Exact BL99 fixed-salinity pieces: FixedDrainedIceSalinityProfile, FixedSalinityBrinePocketEnergyRelation, MaykutUntersteinerConductivity, and BubblyBrineConductivity.
• Descriptive BL99 boundary aliases: MeltingConstrainedSurfaceFluxBalance and OceanFreezingTemperatureBoundary.
• Semi-implicit column energy and salinity stepping, shortwave absorption, budgets, Stefan thickness helpers, melting-limited surface flux handling, and conservative thickness remapping.
• A MutableVerticalDiscretization path that advances energy and prognostic salinity with the moving vertical Jacobian and explicit swept-face enthalpy/salinity fluxes.
• Local CICE/Icepack BL99 column setup, parser, analyzer, generated validation CSVs, and both aggregate and dynamic visualization assets.
• Physics and validation docs for the fixed-salinity relation, moving-coordinate equations, CICE source matching, and current diagnostic limitations.

Moving Coordinate / Changing Thickness

Changing ice thickness is handled with the conservative moving-coordinate balance on Oceananigans' mutable vertical grid. With reference coordinate r, physical coordinate z(r, t), and vertical Jacobian J = partial z / partial r, the implemented energy and salinity equations are

$$\partial_t(J E) = \partial_r(F^E + I) + \partial_r(\dot z_g E^{up}), \qquad \partial_t(J S) = \partial_r F^S + \partial_r(\dot z_g S^{up}).$$

The finite-volume energy update is

$$J^{n+1} \Delta r_k E_k^{n+1} = J^n \Delta r_k E_k^n + \delta z_{g,k+1/2} E^{up}_{k+1/2} - \delta z_{g,k-1/2} E^{up}_{k-1/2} + \Delta t \left[(F^E + I)_{k+1/2}^{n+1} - (F^E + I)_{k-1/2}^{n+1}\right].$$

The assembled tridiagonal solve divides this balance by the current layer thickness, so the right-hand side contains the previous-to-current metric ratio plus the explicit swept-face tracer integral. The salinity step uses the same moving-Jacobian/swept-face structure when bulk salinity is prognostic.

For strict CICE/Icepack BL99 agreement, the validation also applies Icepack's split thickness_changes enthalpy repartitioning for top ablation and basal congelation. The continuous moving-metric replay is still reported as a diagnostic because CICE does not implement growth and ablation as a single continuous moving-grid diffusion solve.

CICE/Icepack BL99 Validation

The validation directory runs and analyzes no-snow one-column CICE histories for:

• fixed cold conductive relaxation
• controlled surface warming
• forced surface ablation
• basal growth

Each case is run with both conduct = 'MU71' and conduct = 'bubbly', for eight histories total.

The aggregate CSV is

validation/cice_bitz_lipscomb/results/cice_case_comparison_summary.csv

All eight histories now have strict_bl99_validation_status = pass. The strict gate uses the source-level BL99 temperature matrix, source-level prognostic-thickness replay, CICE-thickness temperature acceptance, CICE-negative-enthalpy error, and column-energy error. The general ClimaSeaIce split solve and moving-metric CICE-thickness replays are retained as diagnostics; the general path still does not fully reproduce forced top ablation without the source-level Icepack split repartitioning sequence.

Aggregate one-figure summary:

Dynamic thickness-changing visualization for the default basal_growth / MU71 case:

Movie output:

validation/cice_bitz_lipscomb/results/basal_growth_MU71_column_dynamics.mp4

Example Results

The literated example still compares fixed-salinity BL-style and evolving-salinity mushy-column configurations. For the 12-day one-meter comparison with identical forcing and initial profiles:

Column comparison after 12.0 days
  BL surface temperature:    -17.864 C
  Mushy surface temperature: -17.894 C
  Max |T_BL - T_mushy|:      3.071e-02 K
  Max |S_BL - S_mushy|:      1.295e+00 psu
  BL energy budget residual: 2.536e-15
  Mushy energy residual:     0.000e+00
  Mushy salt residual:       1.004e-13

For the melting-limited surface-flux case:

Surface melting case after 2.0 days
  Final thickness:          0.636 m
  Thickness lost:           0.364 m
  Final column temperature: -0.162 C
  Maximum melt flux:        1000.000 W m^-2
  Max budget residual:      0.000e+00

Validation Commands

Ran locally:

• TEST_GROUP=column_energy JULIA_DEPOT_PATH=/private/tmp/julia_depot:/Users/gregorywagner/.julia /Users/gregorywagner/.julia/juliaup/julia-1.11.9+0.aarch64.apple.darwin14/bin/julia --project=. test/runtests.jl
  • column-energy test group passed, including moving-grid metric tests and BL99 public alias checks.
• Regenerated all eight CICE analyzer summaries with validation/cice_bitz_lipscomb/analyze_cice_comparison.jl.
  • all eight aggregate rows report strict_bl99_validation_status = pass.
• JULIA_DEPOT_PATH=/private/tmp/julia_depot:/Users/gregorywagner/.julia /Users/gregorywagner/.julia/juliaup/julia-1.11.9+0.aarch64.apple.darwin14/bin/julia --project=docs validation/cice_bitz_lipscomb/plot_cice_single_column_validation.jl
  • generated validation/cice_bitz_lipscomb/results/cice_single_column_validation.png.
• JULIA_DEPOT_PATH=/private/tmp/julia_depot:/Users/gregorywagner/.julia /Users/gregorywagner/.julia/juliaup/julia-1.11.9+0.aarch64.apple.darwin14/bin/julia --project=. validation/cice_bitz_lipscomb/visualize_cice_column_dynamics.jl
  • generated validation/cice_bitz_lipscomb/results/basal_growth_MU71_column_dynamics.png and validation/cice_bitz_lipscomb/results/basal_growth_MU71_column_dynamics.mp4.
• JULIA_DEPOT_PATH=/private/tmp/julia_depot:/Users/gregorywagner/.julia /Users/gregorywagner/.julia/juliaup/julia-1.11.9+0.aarch64.apple.darwin14/bin/julia --project=. -e 'using ClimaSeaIce; @assert MeltingConstrainedSurfaceFluxBalance === MeltingConstrainedFluxBalance; @assert OceanFreezingTemperatureBoundary(; salinity=34).salinity == 34'
• bash -n validation/cice_bitz_lipscomb/run_cice_cases.sh
• git diff --check

Docs were not rebuilt with docs/make.jl after the latest visualization-doc edit because that script deploys documentation as part of its normal flow.

[glwagner]

Here is a pdf with notes about this model:

A_hierarchy_of_thermodynamic_sea_ice_models.pdf

[glwagner]

@simone-silvestri this PR implements a multilayer sea ice model, following the paradigm proposed in the above PDF. The idea is to solve an energy equation of a particular form that admits a tridiagonal solve for the energy transport. The formulation is general, still, so that existing sea ice models may be implemented into it. The goal here is to show that we can reproduce a few different column implementations.

[simone-silvestri]

Nice! Tomorrow I will go through it!