PaperlayerView exact revision on arXiv

Source-linked AI summary

Physics-enriched neural solvers for transient ice-flow simulation

Thomas Gregov, Sebastian Rosier, Brandon Finley, Andreas Vieli, Guillaume Jouvet

arXiv:2609.12900v1physics.comp-phcs.LGphysics.geo-ph

TL;DR

Transient glacier simulation repeatedly requires solving a nonlinear ice-flow problem as geometry evolves. This paper enriches an online neural solver with inexpensive fields derived from low-order ice-flow balances while leaving its energy objective unchanged. Across three glaciers, the enriched solver is more robust to solver settings and reduces surface-velocity errors by factors of two to four at fixed runtime on the alpine cases.

  • Problem

    Transient glacier simulations repeatedly solve a nonlinear ice-flow problem as glacier geometry evolves.

    Source: §1.1 · p. 1

  • Method

    The solver receives inexpensive fields derived from low-order ice-flow balances as additional network inputs, while the energy, discretization, and boundary conditions remain unchanged.

    Source: §1.3 · p. 3, §1.2 · p. 3, §2.2 · p. 4

  • Results

    Surface-velocity errors fall by factors of two to four at fixed runtime on the alpine cases, while physics-enriched inputs improve solver robustness across all three glaciers.

    Source: §5 · p. 8, §5.1 · p. 9, §6.1 · p. 11

  • Takeaways & Limitations

    Physics-derived inputs better-condition the online minimization, so a few per-solve iterations go further without changing the energy objective.

    Source: §1.2 · p. 3

    The study covers three glacier configurations using Adam, and whether the same rankings hold under quasi-Newton or natural-gradient updates remains open.

    Source: §6.4 · p. 13

Abstract

from arXiv · show

Transient glacier simulations with higher-order ice flow require the repeated solution of a nonlinear problem as the geometry evolves. In the online mode of the Instructed Glacier Model, the velocity field is represented by a neural network whose weights are warm-started from the previous time step and updated with a few optimizer iterations. We show that supplying the network with inexpensive input fields derived from low-order ice-flow balances improves this online solver. Unlike residual-based physics-informed neural networks, which incorporate physics through governing-equation penalties in the loss, our approach leaves the governing energy objective unchanged, adding physical structure through the network inputs. Across three real-world glacier configurations, the enriched solver is markedly more robust to solver settings. On the two alpine cases, it also improves the tuned accuracy--runtime trade-off, reducing surface-velocity errors by factors of two to four at fixed runtime and reaching few-percent relative errors with only $10^4$--$10^5$ trainable parameters, far fewer than comparable raw-input baselines. A 300-year Aletsch simulation then completes in under one minute, and the larger Valais domain in about two minutes, on a single GPU---a budget once reserved for much simpler shallow-ice models. Gains are smaller for the fast marine-terminating glacier, where nonlocal stress coupling favors larger or spectral networks. More broadly, the results suggest that enriching a neural solver's inputs with reduced-order physics can make repeated higher-order solves much cheaper, with no training data and no offline training.

1 Introduction

Higher-order ice-flow simulations repeatedly solve an expensive nonlinear problem, and the online neural solver must make progress with only a few iterations per time step. This paper tests whether inexpensive physics-derived network inputs can make that solver more reliable and efficient without changing its physical objective.

1.1 Solving ice flow with a neural network, online

  • Higher-order ice flow is an expensive nonlinear problem that must be solved repeatedly for glacier predictions, while the online solver warm-starts network weights and applies only a few optimizer iterations per time step.Unlike offline emulators, this solver re-optimizes against the current glacier state's physical energy and uses no training dataset or offline training phase.Source: §1.1 · pp. 1–2

1.2 Physical structure via network inputs

  • DahuNet adds inexpensive, deterministic fields derived from low-order ice-flow balances to network inputs, leaving the energy, discretization, and boundary conditions unchanged.This input-based route adds physical structure where penalty terms in the loss are unavailable or would change the problem being solved.Source: §1.2 · pp. 2–3

1.3 Contributions and outline

  • Across three real-world glacier configurations, physics-enriched inputs improve robustness to solver settings and reduce tuned surface-velocity errors by factors of two to four on the alpine cases.The comparisons include slow alpine Aletsch, the larger Valais domain, and fast marine-terminating Sermeq Kujalleq.Source: §1.3 · p. 3
  • The physics-enriched family is most reliable across the three glaciers, with larger gains on alpine cases; for the marine-terminating case, useful features must match the dominant flow physics.The paper proposes feeding inexpensive limiting balances to networks as inputs when repeatedly solving high-fidelity nonlinear problems along changing trajectories.Source: §1.3 · p. 3

2 An online neural solver for higher-order ice flow

The online solver represents higher-order ice flow as a neural-network parameterization of a discretized energy minimization. It warm-starts weights between nearby glacier states and performs a small, fixed number of optimization iterations at each step.

2.1 Ice flow as energy minimization

  • The Blatter–Pattyn formulation is a nonlinear elliptic problem whose horizontal velocity minimizes a discretized energy containing viscous dissipation, basal friction, gravitational work, and, at marine calving fronts, a water-pressure contribution.The direct minimization over all velocity degrees of freedom serves as the reference problem for the neural networks.Source: §2.1 · p. 4

2.2 The neural solver

  • IGM replaces direct optimization over velocity degrees of freedom with a neural network that maps an input tensor to velocity, then minimizes the same energy over network weights.Raw-input networks receive the glacier state; physics-enriched networks append derived feature fields, while the energy, discretization, and boundary conditions remain identical.Source: §2.2 · pp. 4–5
  • In transient simulations, the solver carries network weights between nearby glacier states and applies a small, fixed number of optimizer iterations to the current state's energy.The resulting velocity advances the glacier state, which in turn defines the energy for the next weight update.Source: §2.2 · p. 5

3 Physics-enriched inputs

DahuNet adds deterministic, inexpensive fields derived from low-order ice-flow balances to the network input rather than changing the energy objective. The study compares raw-input network backbones and enriched variants, including a scrambled-feature control.

3 Physics-enriched inputs

  • Physics enrichment adds deterministic, inexpensive fields derived from low-order balances to the current glacier state, without replacing the higher-order solve or adding trainable parameters.The features are computed with a few elementwise tensor operations per time step.Source: §3 · p. 5

3.1 The feature fields

  • Slope fields expose local flow direction, while the five-feature map also includes velocity-scale fields based on h^4|∇s|^3 and h^3|∇s|^3, with logarithms compressing their wide dynamic range.The fields use fixed, non-learned normalization; the three-feature map contains slope information only.Source: §3.1 · pp. 5–6

3.2 Networks compared

  • The raw-input baselines are a per-pixel MLP, a CNN with local convolutional coupling, and an FNO with global spectral mixing; DahuNet denotes the CNN backbone with five features.CNN+3f, FNO+3f, and FNO+5f are ablation variants.Source: §3.2 · p. 6
  • A scrambled-feature control keeps DahuNet's CNN backbone and channel count but relocates feature patterns by fixed Fourier-phase randomization, testing whether physical feature placement matters.The scrambling preserves field statistics while separating high slopes and large velocity scales from their corresponding locations in the glacier.Source: §3.2 · p. 6

4 Experimental protocol

The protocol evaluates online solver performance through transient-simulation hyperparameter sweeps and uses a separate transfer test to distinguish weight transfer from online correction. It tests three glacier configurations and measures surface-velocity accuracy against direct-solve references.

4 Experimental protocol

  • The evaluation combines hyperparameter sweeps of end-to-end transient simulations with a transfer test on prescribed states from a stored reference trajectory.The transfer test gives every architecture the same states.Source: §4 · p. 6, §4.3 · p. 7

4.1 Glacier test cases

  • The three cases span a slow alpine glacier, a large multi-glacier alpine domain, and a fast marine-terminating glacier; they are not a statistically exhaustive benchmark.Figure 2 maps surface-speed magnitude at maximum reference-trajectory extent, with streamlines, margins, terrain, and the ocean shown.Source: §4.1 · p. 6, Figure 2 · p. 7

4.2 Hyperparameter sweeps

  • Uniform sweeps assess reliability without tuning, while NSGA-II Pareto sweeps approximate the best error–runtime trade-off under a prescribed compute budget.Each trial runs a full transient simulation and records surface-velocity error and wall-clock time per online solve.Source: §4.2 · p. 7

4.3 The transfer test

  • The transfer test separates whether weights remain close to a minimizer across states from how effectively online iterations reduce the energy.Its factorization Go = r Gt separates transfer and online-correction contributions; smaller values of both indicate a good solver.Source: §4.3 · pp. 7–8

4.4 Accuracy metrics

  • Surface-velocity error is the root-mean-square difference from direct-solve reference velocities over the glaciated area, aggregated across evaluation states for the sweep objective.Source: §4.4 · p. 8

5 Numerical results

Physics-enriched inputs improve online-solver robustness and tuned performance across three glaciers, but network size, the source of gains, and useful features vary by flow regime.

5.2 Parameter efficiency

  • On Aletsch, fewer than 5 × 10^3 DahuNet parameters reach about 10% relative error and 4 × 10^4 suffice for 4%, while the plain CNN needs roughly seven times more for comparable accuracy.The MLP never beats the zero-velocity baseline on either alpine glacier, indicating that spatial coupling is indispensable.Source: Figure 4 · p. 9, §5.2 · p. 9

5.1 Robustness and tuned accuracy–runtime performance

  • DahuNet’s median errors were about 23, 12, and 342 m yr−1 on Aletsch, Valais, and Sermeq Kujalleq, 20–50% below the best raw-input baseline on all three.It also makes an untuned solver far less likely to fail.Source: §5.1 · pp. 8–9
  • At fixed runtime, DahuNet reaches 4% relative error (2.3 m yr−1) in about 51 s on Aletsch; on Valais it reaches 6% where raw-input baselines do not, with the best raw-input CNN requiring roughly four times longer.Source: §5.1 · p. 9

5.2 Parameter efficiency

  • Under the no-warm-up protocol, the small-network advantage is absent on Sermeq Kujalleq: the raw CNN is competitive across 10^4–10^5 parameters, and low relative error requires more than about 10^5.The stated explanation is nonlocal membrane-stress coupling, which requires larger receptive fields or globally mixing spectral layers.Source: §5.2 · p. 10

5.3 Attribution of the gains: weight transfer and online correction

  • DahuNet’s online advantage has different origins by regime: better weight transfer on Aletsch, a better from-scratch solution on Valais, and convergence across architectures on the marine trunk.The raw-input FNO solves a single state well from scratch but transfers catastrophically on Aletsch; in the marine case, energy gaps and velocity errors are less tightly coupled.Source: §5.3 · pp. 10–11

5.4 Effect of feature content versus channel count

  • Scrambling feature spatial structure degrades optimization: the control fails from-scratch solves on Aletsch and Sermeq, and its alpine velocity errors are roughly an order of magnitude above DahuNet’s.It also needs an order of magnitude more parameters than DahuNet for comparable Aletsch accuracy.Source: §5.4 · p. 11
  • Useful feature content depends on regime: slope features nearly suffice on Aletsch, five features are preferable on Valais, and the three-feature variant fails while the five-feature variant converges on Sermeq.On Valais, the five-feature map has an energy gap three times smaller and half the velocity error; its added proxy fields supply velocity-scale information needed in the tidewater regime.Source: §5.4 · p. 11

6 Discussion and conclusions

The paper shows that physics-derived inputs improve an online neural ice-flow solver without changing its energy objective, and that the practical gains and useful features depend on glacier regime. It also outlines how small networks make transient higher-order simulations affordable and identifies open questions about optimizer choice.

6.1 Summary of findings

  • By factors of two to four, physics-enriched inputs reduce surface-velocity errors at fixed runtime on the alpine cases and improve robustness to hyperparameter choices across all three glaciers.The experiments change only the network inputs, leaving the energy, discretization, and boundary conditions unchanged.Source: §6.1 · p. 11
  • On Aletsch the enriched inputs mainly improve weight transfer, on Valais they mainly improve the from-scratch solution that transfer preserves, and on the marine case they are a convergence prerequisite while backbone choice matters more.The scrambled control performs worse than no added features, indicating that the gains depend on spatial structure rather than extra channels or parameters.Source: §6.1 · p. 12

6.2 Higher-order ice flow becomes affordable online

  • A DahuNet with a few times 10^4 parameters completes the 300-year Aletsch simulation in under a minute and the larger Valais domain in about two minutes on one consumer GPU, at a few percent relative error.The solver uses about ten Adam iterations per time step; the paper says the small-network regime depends on the physical input features.Source: §6.2 · p. 12

6.3 Guidance for building physics-enriched inputs

  • Physics-enriched inputs can generalize beyond glaciology when inexpensive limiting balances expose solution scales, and the feature content should match the dominant balance without prescribing the solution.The paper recommends small networks whose minimizers are reachable within the per-step iteration budget, fixed a priori nondimensionalization, log compression for wide dynamic ranges, and avoiding divisions by fields that may vanish.Source: §6.3 · p. 12

6.4 Limitations and outlook

  • The study covers three glacier configurations and uses Adam; whether the solver rankings persist with quasi-Newton or natural-gradient updates remains open.The paper notes that optimizer choice should be considered alongside time integration because implicit schemes may permit larger time steps.Source: §6.4 · p. 13

A Physical motivation of the feature fields

The appendix motivates two velocity-scale fields from limiting ice-flow balances: shallow-ice deformation and sliding or membrane-dominated flow. Their geometric factors guide the five-feature map.

A Physical motivation of the feature fields

A.1 Shallow-ice deformation scale

  • For Glen exponent n = 3, shallow-ice deformation motivates the geometric factor h^4|∇s|^3 for busia.The scale follows from combining the depth-dependent shear-stress scale with Glen’s law and integrating the shear rate through the ice column; the constants and Arrhenius factor are not duplicated because A is already in the raw state.Source: §A.1 · p. 13

A.2 Sliding and membrane-flow scale

  • For n = 3, both sliding and membrane-dominated balances motivate the local factor h^3|∇s|^3 for buw.A complete shelf velocity scale would additionally involve a horizontal integration length and boundary conditions.Source: §A.2 · p. 13

A.3 Normalization and robustness choices

The appendix records normalization choices for the proxy fields and identifies the tables containing glacier settings and sweep ranges.

A.3 Normalization and robustness choices

  • The two velocity proxies are log-compressed because their dynamic range would otherwise span several orders of magnitude, while all feature fields use fixed, non-learned scales.The stated scales are k0 = 10^-1 for slopes and h0 = 200 m for thickness.Source: §A.3 · p. 13
  • The driving-stress proxy avoids division by basal friction because that field may vanish and make the division numerically fragile.Instead, buw indicates where driving stress and thickness make large sliding or membrane-dominated velocities plausible; friction remains separately available in the raw input.Source: §A.3 · p. 13
  • Table 2 records glacier-specific settings used in the sweeps and prescribed-state transfer test, while Table 3 lists the hyperparameter ranges used in both sweep types.The proxies do not prescribe velocity; the model’s velocity remains the minimizer of the Blatter–Pattyn energy.Source: §A.3 · p. 13, Table 2 · p. 14, Table 3 · p. 14

B Experimental details

The appendix documents the experimental hardware, glacier configurations, sweep protocols, and shared implementation settings used to reproduce the experiments.

B Experimental details

  • All experimental runs were performed on NVIDIA GeForce RTX 4090 GPUs.The appendix presents these settings as details needed to reproduce the Section 4 experiments.Source: B Experimental details · p. 14

B.1 Glacier configurations

  • Table 2 summarizes the three glacier configurations and their reference trajectories.Source: §B.1 · p. 14

B.2 Sweep protocols and search space

  • Uniform sweeps use 50 trials per architecture and glacier, while Pareto sweeps use up to 500 trials and NSGA-II through Optuna to optimize surface-velocity RMSE and measured wall-clock time per solve.The Pareto search uses population size 20.Source: §B.2 · p. 14
  • The raw-input baselines, DahuNet variants, and scrambled control share the same hyperparameter ranges and use Adam with the stated initialization, activations, and ice-flow discretization settings.The runs use β1 = 0.8, β2 = 0.9995, Glorot-uniform initialization, GELU activations, two vertical velocity layers, bilinear horizontal basis functions, and single-precision arithmetic.Source: §B.2 · p. 14

B.3 Transfer-test settings

Transfer-test configurations were selected under glacier-specific online-time limits, with optimizer budgets held consistent across architectures. The reference energy was computed independently of any network, and Table 4 reports transition medians for the transfer metrics.

B.3 Transfer-test settings

  • Aletsch, Valais, and Sermeq Kujalleq trials were selected under online-time limits of 30 s, 80 s, and 60 s, respectively.For each architecture and glacier, the chosen trial was the lowest-RMSE configuration meeting its limit.Source: §B.3 · p. 14
  • From-scratch solves used 2×104 Adam iterations, while online correction used 10 iterations at a shared learning rate of 5 × 10−5.From-scratch learning rates were glacier-specific but architecture-independent: 2 × 10−4 for Aletsch, 2 × 10−3 for Valais, and 1 × 10−3 for Sermeq Kujalleq.Source: §B.3 · pp. 14–15
  • Table 4 reports medians for from-scratch, transfer, and post-online-iteration gaps, reduction factor, and surface-velocity RMSE.Each metric column is a median over transitions taken independently, so the tabulated median row need not satisfy the per-transition factorization exactly.Source: Table 4 · p. 15
  • The network-independent reference energy used 20,000 Adam iterations followed by at most 5000 L-BFGS iterations, with relative stopping tolerance 10−8.Marginally negative energy gaps can occur because the numerical reference has finite tolerance.Source: §B.3 · p. 15

C Complete transfer-test results and variant sweeps

Across the uniform sweeps, physics-enriched variants outperform the scrambled-feature control, while the explanations for small-network performance differ by flow regime. Local features suit alpine flow, whereas marine flow requires learning nonlocal stress propagation.

C Complete transfer-test results and variant sweeps

  • On Aletsch, FNO models solve from scratch well but have a large transfer gap; on Sermeq Kujalleq, CNN+3f and scrambled variants already fail from scratch.The former points to weak weight transfer, while the latter indicates an inadequate parametrization for the marine regime rather than a transfer failure.Source: C Complete transfer-test results and… · p. 15
  • All physics-enriched variants outperform the scrambled control on every glacier, and CNN+5f has the lowest median error on Sermeq Kujalleq.Figure 6 reports median surface-velocity RMSE per glacier and variant, with interquartile error bars.Source: C Complete transfer-test results and… · pp. 15–16, Figure 6 · p. 16

D Flow locality and the small-network regime

  • On Sermeq Kujalleq under no-warm-up, raw CNN is competitive with DahuNet across 10^4–10^5 parameters, but reaches the low-error front only beyond approximately 10^5.The appendix attributes the contrast to flow locality and identifies optimization budget, rather than representational capacity, as the marine obstacle.Source: D Flow locality and the small-network… · p. 16

D.1 Local versus nonlocal flow

  • At Sermeq Kujalleq, CNN receptive fields must span the coupling length, while FNOs capture nonlocality in one spectral layer but have a parameter floor.DahuNet provides local physical content but still needs iterations to learn the propagation operator.Source: §D.1 · p. 16

D.2 Terminus caveat

Near the calving front, the local driving-stress proxy misses contributions that sustain trunk speed. A fixed pretraining budget extends DahuNet’s low-error parameter range on Sermeq Kujalleq, supporting optimization budget—not representational capacity—as the no-warm-up bottleneck.

D.2 Terminus caveat

  • Near flotation, the local proxy collapses and underestimates the fastest trunk flow; proposed feature treatments still cannot supply the nonlocal grounding-line contribution.The limitation matters more for configurations with extensive floating ice.Source: §D.2 · pp. 16–17

D.3 Optimization budget versus representational capacity

  • With pretraining, the Sermeq Kujalleq DahuNet front extends to a few thousand parameters and is the best option across 10^4–10^5 parameters.Figure 7 grants each trial 5000 pretraining iterations on the first stored state before the transient run.Source: §D.3 · p. 17, Figure 7 · p. 17
Loading 2609.12900v1…