EVP solver fails to converge: ice velocities reach 1.5+ m/s and never return to physical range

[ruoyan88]

Summary

Running a coupled ocean–sea-ice simulation on a 0.25° Arctic regional grid (55°N–90°N, TripolarGrid), sea-ice velocities grow monotonically from rest to ~1.5 m/s within the first 6 hours and remain at that unphysical level indefinitely. Halving the time step from 2 min to 1 min has no meaningful effect on ice velocity magnitudes, suggesting the EVP rheology is not generating sufficient internal stress to arrest ice motion. Typical Arctic sea-ice drift speeds are 0.05–0.3 m/s.

Environment

Component	Version
Julia	1.10.10
Oceananigans	0.104.2
ClimaSeaIce	0.4.2
ClimaOcean	0.9.0
CUDA.jl	5.9.6
CUDA runtime	12.9
CUDA driver	550.163.1 (supports up to 12.4)
GPU	4 × NVIDIA GeForce RTX 4090 (24 GB each)

Halving Δt improves thickness conservation but has essentially zero effect on the velocity magnitude, confirming this is not a CFL violation but an EVP convergence failure.

Expected Behavior

With 500 EVP substeps and Δt = 60 s (Δt_evp = 0.12 s), the elastic waves should damp sufficiently for the stress state to converge toward the VP solution. Ice velocities should stabilize at O(0.1) m/s within the first few hours of integration, consistent with the applied atmospheric forcing and realistic ice rheology.

Hypotheses

1. Insufficient EVP damping: 500 substeps may not be enough for this grid resolution (Δx ≈ 7–14 km at 60°–80°N). The elastic wave CFL for EVP may require more substeps or different damping parameters (e.g., mEVP α/β tuning).

2. Ice pressure too low at initialization: If the VP ice pressure P = P* · h · exp(−C(1−ℵ)) is not properly activated during the first few steps (e.g., because strain rates are initially zero), the ice enters free-drift and builds up too much kinetic energy before internal stress can respond.

3. Spinup shock with strong atmospheric forcing: Starting from rest with full JRA55 wind forcing may overwhelm the EVP solver before it has time to build up internal stress. A gradual ramp-up of atmospheric forcing might help, but this should not be required for a stable solver.

4. Possible issue with stress boundary conditions at sponge edges: The velocity maximum at exactly φ = 60°N (sponge boundary) suggests the restoring force may interact adversely with EVP stress computation.

Questions for developers

1. What EVP variant does ClimaSeaIce use (standard EVP, mEVP, or aEVP)? What are the recommended substep counts and damping parameters for 0.25° resolution?

2. Is there a known issue with EVP convergence when starting from rest with strong forcing?

3. Are there any internal velocity limiters or adaptive damping mechanisms we should enable?

4. For the elastic parameters (E or T_evp), what values does the code use by default, and are they appropriate for Δt = 60 s with 500 substeps?

Reproduction

Happy to provide the full simulation script and initial condition files. The issue is reproducible on every run from the same IC with both Δt = 60 s and Δt = 120 s. All structural verification checks (grid, sponge coverage, array dimensions, broadcasting) pass cleanly.

Log excerpt

First 12 hours of STEP 5 velocity monitoring

iter:   0, time: 0 s,       Δt: 1 min | uᵢ: (0.00, 0.00) m/s | uₒ: (0.00, 0.00) m/s | h: 4.12 m
iter:  10, time: 10 min,    Δt: 1 min | uᵢ: (0.46, 0.41) m/s | uₒ: (0.11, 0.10) m/s | h: 4.12 m
iter:  40, time: 40 min,    Δt: 1 min | uᵢ: (0.55, 0.60) m/s | uₒ: (0.43, 0.41) m/s | h: 4.12 m
  ⚠ [STEP5] iter=40: ice velocity continuously increasing (0.55 m/s), check EVP convergence
iter:  60, time: 1 h,       Δt: 1 min | uᵢ: (0.58, 0.69) m/s | uₒ: (0.62, 0.60) m/s | h: 4.12 m
iter: 120, time: 2 h,       Δt: 1 min | uᵢ: (0.77, 0.93) m/s | uₒ: (0.98, 1.05) m/s | h: 4.12 m
iter: 220, time: 3.67 h,    Δt: 1 min | uᵢ: (1.01, 1.11) m/s | uₒ: (1.35, 1.28) m/s | h: 4.11 m
  ⚠ Ice velocity exceeding 1 m/s
iter: 240, time: 4 h,       Δt: 1 min | uᵢ: (1.17, 1.15) m/s | uₒ: (1.42, 1.33) m/s | h: 4.11 m
iter: 360, time: 6 h,       Δt: 1 min | uᵢ: (1.55, 1.55) m/s | uₒ: (1.66, 1.32) m/s | h: 4.10 m
iter: 480, time: 8 h,       Δt: 1 min | uᵢ: (1.50, 1.45) m/s | uₒ: (1.45, 1.36) m/s | h: 4.09 m
iter: 600, time: 10 h,      Δt: 1 min | uᵢ: (1.49, 1.54) m/s | uₒ: (1.75, 1.31) m/s | h: 4.09 m
iter: 700, time: 11.67 h,   Δt: 1 min | uᵢ: (1.48, 1.61) m/s | uₒ: (1.70, 1.44) m/s | h: 4.09 m

[simone-silvestri]

Hello @ruoyan88. To answer your questions:

1. ClimaSeaIce uses an aEVP formulation.
2. I don't think so, I did not see any unphysical behavior happening
3. No they are not.

What I wanted to ask is if you are sure the velocities were unphysical. Note that ClimaSeaIce computes the velocities also where there is no sea ice, which fall back to a balance of stresses (basically to the ocean surface velocities). This might be where you have u ~ 1.5 m / s. You probably want to check u * \aleph, where \aleph is the ice concentration.

[glwagner]

@ruoyan88 can you please include a plot or animation of your solution? Also, we will need your entire script to reproduce the issue. Without a script that reproduces the issue, there is little action we can take.