Leveraging higher-order time integration methods for improved computational efficiency in a rainshaft model

Abstract

Cloud and precipitation microphysics packages in atmospheric general circulation models typically use first-order time integration methods with a large time step, requiring ad hoc limiters and substepping of the sedimentation scheme to prevent solutions from becoming unstable. We show that in the latest version of Energy Exascale Earth System Model, E3SMv3, the rain microphysics provided by the Predicted Particle Properties (P3) scheme is underresolved in time at the model's default 300s time step. The P3 scheme requires limiters to guarantee stability, but those limiters make large discretization errors more difficult to detect. When the time step of the P3 scheme is reduced to sufficiently capture correct microphysics behavior, wall clock time of the simulation is increased by nearly a factor of 40. Instead of reducing the microphysics time step, we recommend using higher-order time integrators based on Runge-Kutta methods, which offer improved solution accuracy at comparable computational costs. A key to obtaining computationally efficient microphysics results is the use of adaptive time stepping, which also eliminates the need for specialized substepping procedures in the sedimentation process. We also analyze individual microphysical processes by extracting inverse timescales from Jacobians of the process rates, which gives insight about the maximum time step each process is able to take while maintaining stability and accuracy, and about how individual processes should be grouped together for most efficient results. The proposed integrators can achieve the accuracy level required to correctly model rain microphysics parameterizations more than 10x faster than the P3 scheme.

Figures (13)

• Figure 1: 2D histogram of pointwise log difference between highly-resolved reference and P3 solution with time steps of 300s, 3.2768s, and 0.4096s, with the differences shown as percentages relative to the highly-resolved reference solution.
• Figure 2: Selected solution profiles at 300s using the P3 time integration method with $\Delta t = 300s$, 3.2768s, and 0.4096s.
• Figure 3: Left: initial solution profiles for temperature, water vapor mixing ratio, raindrop number concentration, and rain mass mixing ratio obtained from a single column of input to P3. Right: corresponding fall speed $v_{q_{\text{r}}}$ and characteristic speeds $\sigma_{\text{sed},1}$ and $\sigma_{\text{sed},2}$. The fall speed used in P3's substepping calculations $(v_{q_{\text{r}}})$ is notably slower than the speed dictated by the largest sedimentation eigenvalue.
• Figure 4: Inverse timescales for sedimentation, self-collection, and evaporation for 1000 columns from E3SM described in Section \ref{['sec:boundary conditions']} at $t=0s$ and $t=200s$. Generally, sedimentation timescales are faster than those of the other processes, although in some columns self-collection timescales are as fast as or faster than those of sedimentation. Additionally, the sedimentation timescale for stability ($\sigma_{\text{sed},2}/\sqrt{\frac{d^{2}q_{\text{r}}}{dz^{2}} / q_{\text{r}}}$) is faster than the sedimentation timescale for accuracy ($\sigma_{\text{sed},2}/\Delta z$), suggesting that sedimentation is a mildly stiff process.
• Figure 5: Predicted stability bound on time step for column 37 based on (i) CFL condition (see Equation \ref{['eq:cfl']}), (ii) fall speeds (see Equation \ref{['eq:P3 cfl']}), and (iii) stability region of the forward Euler method when considering the rainshaft model (Equation \ref{['eq:rainshaft model']}) with sedimentation as the only active process. We compare these bounds with the time step selected by an adaptive time step approach in SUNDIALS with a second-order ERK method which uses an embedded first-order method to estimate temporal error and adapt the time step. The SUNDIALS controller performs well in correctly capturing the CFL condition (Equation \ref{['eq:cfl']}).