Neural ODEs and symbolic regression recover an ice-crystal growth equation from 290 mass time series and reproduce early growth in independent AIDA chamber data
Synopsis
Using neural ordinary differential equations (NODEs) optimized simultaneously across 290 mass-ratio time series of ice crystals grown in a levitation diffusion chamber, the study learned the unknown functional form of the transfer coefficient G in the depositional ice growth model, then used symbolic regression (PySR) to derive a closed-form expression; the weakly-constrained NODE model performed best on 138 of 290 experiments (MSE loss 16495, versus 41396 for no surface kinetics and 40751 for Nelson and Baker 1996), yielding G = a0·Gc^a1·[a2 + a3/m]^-1 + a4 with an added ice-crystal-mass term, and this expression predicted ice supersaturation and ice water content better than the capacitance growth model and Nelson and Baker 1996 in independent AIDA Aerosol and Cloud Chamber experiments.
Figure 1: Overview of methodology for learning unknown physics in depositional ice growth models. Depositional ice growth is an important process for ice formation in atmospheric clouds. Ice crystals grow via direct deposition of water molecules from the vapor phase onto the ice surface. a) We replace partially unknown physics in the depositional ice growth model with a neural network, considering both a strong and weak constraint. b) We integrate the ice growth rate partially parameterized by a neural network and optimize to reduce the distance between the model and experimentally measured time series of ice mass ratios. c) We use symbolic regression to determine a functional form for the unknown physics learned by the neural network.
· Page 3Interpretation
The weakly-constrained NODE model (treating G/Gc as a neural network function of Si, T, and m) achieved the lowest MSE loss on 138 of 290 levitation diffusion chamber experiments, with an overall MSE loss of 16495, outperforming the no-surface-kinetics model (41396), Nelson and Baker 1996 (40751), and the strongly-constrained NODE model (30714). Prior depositional ice growth models relied on parameterizing a single deposition coefficient α, and α values from different experiments disagreed by orders of magnitude; this work instead directly learns the ratio of the transfer coefficient G to the continuum-limit Gc and optimizes across all 290 time series simultaneously, bypassing prior assumptions about the α functional form. Based on 290 1 Hz mass-ratio time series, with the first 500 seconds used for training, compared via MSE loss and per-experiment best-count metrics; the authors also validated the method on a synthetic dataset with a known α functional form, showing it recovers the true nonlinear dependence.
Symbolic regression on the trained neural network yielded the closed-form expression G = a0·Gc^a1·[a2 + a3/m]^-1 + a4 (a0=688.267, a1=1.3153, a2=0.85601, a3=2.6606×10^-12, a4=0.1123×10^-9); all candidate expressions depend on Gc and ice crystal mass m, and only the most complex expressions show any dependence on Si and T. This expression adds a term related to ice crystal mass m to the capacitance growth model, whereas classical capacitance theory describes growth only through Gc; the authors note this is consistent with physical expectations that surface kinetic effects suppress growth for small crystals, and that most of the variance can be attributed to ice crystal size rather than Si or T. PySR was run for 1000 generations using m, r, T, Si, and Gc as input features to fit G, producing a Pareto front of candidate expressions (13 in Table S2); the expression balancing complexity and generalization was selected using independent AIDA data.
On independent data from AIDA Aerosol and Cloud Chamber IsoCloud experiments (195–235 K, 150–300 hPa, average ice crystal sizes <10 µm), the new model predicted ice supersaturation and ice water content better than the capacitance growth model and Nelson and Baker 1996, particularly at lower temperatures. Prior models were mostly validated under single controlled conditions; this work directly substitutes the formula learned from the levitation diffusion chamber into the depositional growth term of an AIDA bin microphysics model, demonstrating transferability to time-evolving realistic cirrus conditions. AIDA experiments are an independent dataset whose Si and T range partially overlaps with the levitation diffusion chamber but includes lower T and higher Si with time-evolving environmental conditions; comparison used a bin microphysics model constrained to observed ice number concentrations and evaluated MSE loss.
The strongly-constrained NODE model (assuming a single deposition coefficient α = fα(Si, T)) failed to learn a consistent α function for most experiments, with predicted α differing significantly from Nelson and Baker 1996 and per-experiment performance only slightly better than the no-surface-kinetics model. This suggests that representing surface kinetic effects as a single deposition coefficient α function may be too restrictive, consistent with Pokrifka et al.'s finding that newly formed crystals require transformations in surface growth modes. Based on MSE loss across 290 experiments (30714) and per-experiment best-count (61, close to the no-surface-kinetics model's 62), plus visualizations of mass-ratio deviations and α distributions in Supplementary Figures S5 and S6.
Perspective
The results target parameterization development for early ice crystal growth in cirrus and upper-level clouds, applicable to the 205–240 K and ice supersaturation 1.0–1.8 conditions covered by the levitation diffusion chamber, and to the 195–235 K, 150–300 hPa cirrus simulation conditions of the AIDA experiments. For climate model developers, Eq. 6 can directly replace the G term in the capacitance growth model; for laboratory researchers, the NODE-plus-symbolic-regression workflow can be applied to parameterization derivation for other microphysical process rates. The authors also note an alternative approach is to derive a distribution of growth rates from these experimental measurements rather than a single functional form.
The authors note that the extrapolation test beyond 500 seconds is biased toward compact, less efficiently growing ice crystals (typically at lower supersaturations), and short time series from rapid polycrystalline growth at high supersaturation may be underrepresented; when integrating well beyond 1000 seconds, the weakly-constrained NODE model begins to over-predict growth at larger sizes, suggesting depositional growth should asymptote to the continuum limit at large mass. Additionally, mass derivation in the levitation diffusion chamber assumes ice crystals are initially spherical, whereas they may be polycrystalline at nucleation; supersaturation uncertainty is about 10% and may show a one-directional bias. These are scope considerations for readers applying the formula to broader conditions.