| 21 Jan 2026
Ocean–atmosphere turbulent flux algorithms in Earth system models do not always converge to unique and physical solutions: analysis and potential remedy in E3SMv2
Abstract. The development of physics parameterizations in Earth system models typically emphasizes whether the intended physics is reasonably represented, while mathematical aspects such as solvability of the governing equations and convergence of the numerical algorithms used to approximate their solutions receive far less attention. In this paper, we examine these mathematical issues for a widely used ocean–atmosphere turbulent flux parameterization and its implementation in the Energy Exascale Earth System Model version 2 (E3SMv2). We show that, under simulated meteorological conditions, the parameterization can yield no solution or multiple (including unintended) solutions. These problems arise primarily from (1) a discontinuity in the formulation of the neutral exchange coefficients and (2) the use of an ad hoc limiter on the Monin–Obukhov length to address a singularity in its definition. Compounding these problems is the fact that interventions of calculations such as limiters are often thought to have only a "minor" effect on numerical algorithms and are not documented in technical model descriptions. To address these solvability issues, we propose (1) a regularization that enforces continuity in the neutral exchange coefficients and (2) an adaptive procedure for selecting limiting values of the Monin–Obukhov length based on mathematical analysis of solution uniqueness. Implementing these revisions in E3SMv2 leads to statistically significant changes in the simulated latent heat fluxes over the mid-latitude oceans in the winter hemisphere as well as over the subtropical and tropical oceans. Overall, this work improves the well-posedness and numerical accuracy of ocean–atmosphere turbulent flux calculations in E3SMv2. Moreover, because discontinuities and ad hoc limiters are frequently encountered in physics parameterizations, this work serves as an example of how non-existence and non-uniqueness issues in parameterizations can be identified, analyzed, and resolved.
General Comments
This is a very welcome paper that analyzes the well-posedness of an atmosphere-ocean surface flux parameterization widely used in climate models (Large and Pond). The paper is within the scope of GMD in that it is a technical paper concerned with analyzing and improving the numerical scheme used to solve the equations in the parameterization. As the authors state, this is likely the first such analysis of the Large and Pond scheme or any turbulent flux scheme used in climate models. The paper is well organized and makes a clear case that changes should be made to existing models using this parameterization. Two flaws in the numerical solutions are found: one related to the discontinuity in the formulation of the neutral exchange coefficient and the other related to the use of an ad-hoc limiter on the Monin-Obhukhov length. These flaws can lead to physically unrealistic fluxes from readily encountered meteorological conditions. Fixes are proposed and examined in specific sets of conditions and in global atmosphere-only simulations with E3SMv2. I share the hope with the authors that this work will inspire similar examinations of algorithmic implementations of other well-known parameterizations and cause model developers to think more carefully about the introduction of ad-hoc limiters.
Specific Comments
The presentation and explanation of the problem and solution is very clear. Just enough background is provided in the text for a climate scientist to understand the mathematical analysis.
When instantaneous daily values of inputs for the flux scheme are saved, is that done at the same simulated time of day? Is it 0Z? And for all 10 years? You might miss some meteorological conditions since half the planet is always night and half always day if always sampling at 0Z but I don't think that would change the results. You should mention how many sets of values were saved. (10 * 365 * # of ocean grid points) total. Was every set subject to the offline calculation mentioned in 3.2?
In section 5, results are presented using a 10-year mean calculated from a 10-year simulation. Even in an F-case, there will be some adjustment time because of the land surface. You should look at a time series of global annual average surface air temperature to see how quickly it reaches a quasi-steady state. You may have to discard the first year and use 9-year averages. If a large adjustment to the new scheme is in the first year of SENS, it would distort the tests of significance.
Technical Comments
I did not notice any glaring spelling errors or grammar issues.
In equation A11, should it be CH(U/u*) and CE(U/u*) ? Something isn't right if √CD = u*/U.