Numerical coupling of aerosol emissions, dry removal, and turbulent mixing in the E3SM Atmosphere Model version 1 (EAMv1) – Part 2: A semi-discrete error analysis framework for assessing coupling schemes

Part 1 (Wan et al., 2024) of this study discusses the motivation and empirical evaluation of a revision to the aerosol-related numerical process coupling in the atmosphere component of the Energy Exascale Earth System Model version 1 (EAMv1) to address the previously reported issue of strong sensitivity of the simulated dust aerosol lifetime and dry removal rate to the model's vertical resolution. This paper complements that empirical justification of the revised scheme with a mathematical justification leveraging a semi-discrete analysis framework for assessing the splitting error of process coupling methods. The framework distinguishes the error due to numerical splitting from the error due to the time integration method(s) used for each individual process. Such a distinction results in a framework that provides an intuitive understanding of the causes of the splitting error. The application of this framework to the dust life cycle in EAMv1 confirms (i) that the original EAMv1 scheme artificially strengthens the effect of dry removal processes and (ii) that the revised splitting reduces that artificial strengthening.

While the error analysis framework is presented in the context of the dust life cycle in EAMv1, the framework can be broadly leveraged to evaluate process coupling schemes, both in other physical problems and for any number of processes. This framework will be particularly powerful when the various process implementations support a variety of time integration approaches. Whereas traditional local truncation error approaches require separate consideration of each combination of time integration methods, this framework enables evaluation of coupling schemes independent of particular time integration approaches for each process while still allowing for the incorporation of these specific time integration errors if so desired. The framework also explains how the splitting error terms result from (i) the integration of individual processes in isolation from other processes and (ii) the choices of input state and time step size for the isolated integration of processes. Such a perspective has the potential for the rapid development of alternative coupling approaches that utilize knowledge both about the desired accuracy and about the computational costs of individual processes.

1 Introduction

Accurate representation of process interactions is an important and ubiquitous challenge in multiphysics modeling (Keyes et al., 2013). For the sake of tractability, splitting methods are widely used that allow for the separate development of both the continuum and discrete representation of individual processes, which are then assembled to form a multi-process numerical model. In weather, climate, and Earth system modeling, it has been recognized that how to combine the different process representations to form a coherent and accurate numerical model is a challenge deserving more attention; see the review paper by Gross et al. (2018) and the references therein as well as more recent studies by authors such as Barrett et al. (2019), Donahue and Caldwell (2020), Wan et al. (2021), Santos et al. (2021), Ubbiali et al. (2021), and Zhou and Harris (2022). The dust aerosol life cycle problem discussed in the companion paper (Wan et al., 2024) is a recent example from the atmosphere component of the Energy Exascale Earth System Model version 1 (EAMv1; Rasch et al., 2019) that shows that different numerical methods used for process coupling in the overall time integration can lead to substantially different results at a fixed spatial resolution as well as to significantly different sensitivities to spatial resolution change.

That companion paper (Wan et al., 2024) clarified the main source and sink processes in the dust life cycle in EAMv1, quantified their relative magnitudes, reflected on the process coupling scheme used in the default model, proposed a revised coupling scheme, and evaluated the impact of the revised coupling on the simulated aerosol climatology. The discussions therein are based primarily on the intuition of atmospheric modelers, and the reasoning was verified by confirming agreement between the expected and obtained numerical results from EAMv1. To gain more confidence that the EAMv1 solution obtained with the revised coupling scheme is indeed a better numerical solution (i.e., closer to the true or trusted solution than that of the original EAMv1 scheme), a mathematical explanation for the changes with the revised coupling is needed. Typically, computational analyses are done to develop such explanations, including time step self-convergence studies, such as those in Wan et al. (2013), or time step sensitivity studies, such as those in Wan et al. (2021) and Santos et al. (2021). Both of these approaches are unfortunately impractical in this case, as (in the current EAM code) the coupling time step for dust emissions, dry removal, and turbulent mixing is tied to the coupling time steps of various other atmospheric processes, such as aerosol microphysics and gas-phase chemistry, deep convection and aerosol wet removal, and the coupling between the resolved dynamics and the parameterizations. Thus, without significant code structure changes, it is not feasible to isolate the impact of coupling approaches for the three aerosol processes that we would like to focus on.

This paper provides a theoretical explanation of the numerical results presented in Part 1 (Wan et al., 2024) and introduces a framework based on truncation error analysis that can more broadly help address the research gap in numerical process modeling in weather, climate, and Earth system modeling. In a fully discretized model with process splitting, the overall local truncation error from time integration includes contributions from (i) the time integration of each individual process and (ii) the splitting of each process from the remaining processes. In the atmosphere modeling literature, there are a number of theoretical studies that leverage truncation error analysis to compare the accuracy of different splitting methods (e.g., Caya et al., 1998; Staniforth et al., 2002; Dubal et al., 2004, 2005, 2006; Ubbiali et al., 2021). Considering the added complexity of the weather, climate, and Earth system models, especially the long-term and large-team efforts that are typically needed for the continuous development of such models, it is useful to understand and reduce the different types of error separately. Because the numerical results in Part 1 (Wan et al., 2024) showed substantial sensitivity to the process coupling approach, this work develops a framework that leverages exact time integration of the processes to focus on understanding and reducing the splitting error. Because the time discretization is also an error source worth assessing and addressing, the framework is designed to be flexible enough to incorporate the individual process time integration errors (see discussion at the end of Sect. 2.4).

In the work of Williamson (2013) that discusses issues related to atmospheric convection and process splitting, exact time integration is also used for individual processes in a two-process equation and a three-process equation. However, the purpose there was to point out a problem in the formulation of a specific parameterization; hence, understandably, that work did not use exact time integration as a general tool for analyzing splitting errors in physics systems beyond the two highly idealized and customized equations discussed therein. In this paper, we demonstrate the usage of exact time integration within an error analysis framework that distinguishes splitting error and the error resulting from temporal discretization of individual processes, providing an approach that allows for a focus on the splitting error in general problems involving two or more processes. While such a focus is commonplace in the field of mathematics (e.g., Hairer et al., 2006; LeVeque, 1982), the approach has not caught much attention in the weather and climate modeling communities, despite the important benefits that its adoption can provide to the development of sophisticated atmospheric numerical models.

For example, in EAM and its predecessors, changes in time integration methods have been implemented both in the dynamical core representing the resolved fluid dynamics and in various parameterizations representing sub-grid-scale processes. The results from a local truncation error analysis that did not distinguish between splitting and time integration sources would become invalid as soon as any time integration method was changed; hence, separate considerations would be needed for every combination of time integration methods used by the dynamical core and the many physics parameterizations. In contrast, the results from an error analysis approach that can consider splitting error independent of time integration sources are expected to have a better chance of remaining valid across multiple versions of the same atmosphere model and might even generalize to other atmosphere models.

Another benefit of the splitting error analysis framework demonstrated in this work is that the terms in error expressions produced are easily attributed to the coupling choices made for the individual right-hand-side (RHS) terms of the continuum equations. In atmosphere modeling terminology, the framework goes beyond deriving the splitting errors that contribute to the overall error in prognostic variables; it also demonstrates how the coupling choices lead to errors in the process rates (i.e., rates of change in prognostic variables, also referred to as tendencies) associated with the various physical processes. For this work, coupling choices for the dust source and sink processes are directly mapped to the splitting error terms that contribute to overall error in the EAM-simulated mixing ratios. Reducing splitting errors in the process rates can help avoid compensating errors from different physical processes and, thus, help ensure the model provides good predictions of the prognostic variables for the right reasons. Additionally, understanding the impact of numerical coupling scheme choices at the process level allows for the development of new coupling strategies that focus on improving the accuracy of numerical process rates (and the associated prognostic variables) while considering the computational cost of the various processes.

Whereas the companion paper (Wan et al., 2024) focuses on motivating the dust life cycle problem and an empirical comparison of the two coupling methods in consideration, this paper focuses on how a semi-discrete error analysis supports the empirical finding that one coupling method leads to a better numerical solution than the other. This two-part approach facilitates a detailed description of the framework with significant pedagogical value to weather and climate model developers. For example, our own collaboration between applied mathematicians and atmospheric scientists on the investigation of dust life cycle in EAMv1 has shown that a step-by-step explanation of the derivation of splitting errors resulting from two coupling methods, which are widely used in weather and climate models, was helpful for the EAM developers in this collaboration to recognize the relevance, as well as the generality, of the semi-discrete methodology. In addition, one of the points that we make in Sect. 3 is that, after splitting errors are derived for coupling schemes used in two-process problems, it is possible to use those error expressions as building blocks to perform back-of-the-envelope derivations for problems involving more processes, making the derivations much less tedious and the framework much easier to use by researchers of specific applications. The discussion in Sect. 3.2 can be viewed as an example of such a back-of-the-envelope derivation, demonstrating that the mathematically rigorous framework can be made accessible to physical scientists working on practical problems. Given the reinvigorated interest in numerical process coupling reflected in the review by Gross et al. (2018) and the community efforts described in publications such as Heinzeller et al. (2023), the pedagogical description of the error analysis framework presented here can be a useful contribution to those model development efforts.

To focus on the error analysis framework, the remainder of this paper forgoes the background on the motivating dust life cycle problem and coupling methods in consideration, for which the reader is referred to Part 1 of this work (Wan et al., 2024), and instead opens by introducing the analysis framework in Sect. 2 using a generic two-process problem and deriving the splitting errors associated with the widely used parallel and sequential splitting methods. Section 3 then applies the error analysis framework to a three-process problem inspired by the dust aerosol life cycle in EAMv1. Leading-order splitting errors are derived for both the original process coupling in EAMv1 and the revision proposed in Part 1, and the characteristics of the splitting errors are discussed. This paper concludes by summarizing results in Sect. 4. Appendix A and B provide mathematical details of the analysis framework.

2 A semi-discrete analysis framework for assessing splitting error of process coupling methods

The error analysis framework demonstrated in this work is described as semi-discrete in that it takes the perspective that, while the overall time integration of the model is discrete, the integration of individual processes is done exactly (i.e., there is no temporal discretization for each process). This perspective is the critical component that allows the framework to isolate the splitting truncation error from that of the process temporal discretization errors. The semi-discrete approach allows for the derivation of splitting truncation errors from how the coupling scheme incorporates individual process. By casting the numerical splitting of processes as estimating process rates in a split fashion, the framework identifies two general coupling method choices that cause splitting truncation errors. In Sect. 2.2, those coupling choices are identified in two widely used coupling methods for two-process problems. In Sect. 2.3, the splitting truncation error terms that result from those choices are derived and combined to form the leading-order splitting truncation error for each coupling method. In Sect. 2.4, the terms in the leading-order splitting truncation error are attributed back to the coupling method choices in a manner that (i) identifies when the splitting truncation errors from coupling methods compound or cancel each other and (ii) provides a workflow to easily generalize the results to problems with more than two processes.

2.1 Notation and definitions

Consider a prognostic equation, with multiple operators, that can either be an ordinary differential equation (ODE) or a partial differential equation (PDE). Noting that the “method of lines” will reduce a PDE to an ODE, the prognostic equation to be studied herein is introduced in the ODE form:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}q\left(t\right)}{\mathrm{d}t}}=\sum _{i=\mathrm{1}}^I{X}_i\left(q\right(t\left)\right),I>\mathrm{1},t>\mathrm{0},q\left(\mathrm{0}\right)=q^{\mathrm{IC}},\end{array}$$

where the different processes, X_i, are discretized and implemented by different components of the model software (i.e., the processes are split). At time t_n, denote q(t_n) and qⁿ as the exact and numerical (approximate) solutions, respectively. The error at time t_n is defined as

$${\mathcal{E}}^n\equiv q^n-q\left(t_n\right).$$

Let ℱ_Δt(q^input) represent the numerical algorithm that advances the solution from state q^input to state q^output, where Δt is the time step size used. In other words, denote ℱ_Δt as the mapping such that the numerical solution at $t_{n+\mathrm{1}}=t_n+\mathrm{\Delta }t$ is

$$q^{n+\mathrm{1}}={\mathcal{F}}_{\mathrm{\Delta }t}\left(q^n\right).$$

The solution error at time t_n+1 can be expressed as

$$\begin{array}{ll}{\displaystyle }{\mathcal{E}}^{n+\mathrm{1}}=& {\displaystyle }{\mathcal{F}}_{\mathrm{\Delta }t}\left(q^n\right)-q\left(t_{n+\mathrm{1}}\right)\\ {\displaystyle }=& {\displaystyle }{\mathcal{F}}_{\mathrm{\Delta }t}\left(q^n\right)-{\mathcal{F}}_{\mathrm{\Delta }t}\left(q\right(t_n\left)\right)+{\mathcal{F}}_{\mathrm{\Delta }t}\left(q\right(t_n\left)\right)-q\left(t_{n+\mathrm{1}}\right)\\ \text{(2)}& {\displaystyle }=& {\displaystyle }\underset{\text{propagated error}}{\underbrace{{\displaystyle \frac{\mathrm{d}{\mathcal{F}}_{\mathrm{\Delta }t}}{\mathrm{d}q}}\left({\mathit{\gamma }}_n\right){\mathcal{E}}^n}}+\underset{\text{local truncation error }\left(\mathrm{lte}\right)}{\underbrace{{\mathcal{F}}_{\mathrm{\Delta }t}\left(q\right(t_n\left)\right)-q\left(t_{n+\mathrm{1}}\right)}},\end{array}$$

where γ_n is a value of q between q(t_n) and qⁿ by the mean value theorem. The error ℰⁿ⁺¹ consists of the evolution of the existing error at time t_n (the propagated error) and the generation of new error from time t_n to t_n+1 (the local truncation error).

If each aforementioned component of the model software implements a single process in Eq. (1) without using information about other processes, then the role of such a component of the software can be interpreted as integrating the following one-process ODE:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{d}{q}_{X_i}\left(t\right)}{\mathrm{d}t}}=X_i\left(q_{X_i}\right(t\left)\right),t>t_n,q_{X_i}\left(t_n\right)=q_{X_i}^{\mathrm{input}}.\end{array}$$

Depending on which coupling scheme is used, $q_{X_i}^{\mathrm{input}}$ can be the numerical solution of the multi-process problem at t_n (i.e., qⁿ) or some value of the physical quantity, q, passed to the component of the software, X_i, by another component of the software, X_j. In the following, the notation $q_{X_i}[t-t_n;q_{X_i}^{\mathrm{input}}]$ is used to denote the exact solution of the one-process ODE (Eq. 3), namely,

$$\begin{array}{ll}{\displaystyle }{q}_{X_i}[t-t_n;q_{X_i}^{\mathrm{input}}]& {\displaystyle }\equiv q_{X_i}\left(t\right)\\ \text{(4)}& {\displaystyle }& {\displaystyle }=q_{X_i}^{\mathrm{input}}+\underset{t_n}{\overset{t}{\int }}{X}_i\left(q_{X_i}\left(\mathit{\eta }\right)\right)\mathrm{d}\mathit{\eta }.\end{array}$$

The explicit mentioning of $q_{X_i}^{\mathrm{input}}$ in the square brackets on the left side of Eq. (4) emphasizes the dependence of $q_{X_i}$ and X_i on $q_{X_i}^{\mathrm{input}}$, the significance of which will become clear below. Using a similar notation for the exact solution of the multi-process problem Eq. (1), one can write

$$\begin{array}{ll}{\displaystyle }q[t-t_n;q(t_n\left)\right]& {\displaystyle }\equiv q\left(t\right)\\ \text{(5)}& {\displaystyle }& {\displaystyle }=q\left(t_n\right)+\underset{t_n}{\overset{t}{\int }}\sum _{i=\mathrm{1}}^I{X}_i\left(q\right(\mathit{\eta }\left)\right)\mathrm{d}\mathit{\eta }.\end{array}$$

To facilitate comprehension of the derivations below, it is again emphasized that the time integrals in Eqs. (4) and (5) are assumed to be exactly evaluated. It is also useful to note that, by definition,

$$\begin{array}{}\text{(6)}& {\displaystyle }{q}_{X_i}[\mathrm{0};q_{X_i}^{\mathrm{input}}]=q_{X_i}^{\mathrm{input}},q[\mathrm{0};q(t_n\left)\right]=q\left(t_n\right).\end{array}$$

2.2 Sources of splitting error

From Eq. (5), one can see the average process rate for X_i that contributes to the change in q from t_n to t_n+Δt in the original multi-process ODE is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\underset{t_n}{\overset{t_n+\mathrm{\Delta }t}{\int }}{X}_i\left(q\right(\mathit{\eta }\left)\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\underset{t_n}{\overset{t_n+\mathrm{\Delta }t}{\int }}{X}_i\left(q\right[\mathit{\eta }-t_n;q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}{X}_i\left(q\right[\stackrel{\mathrm{̃}}{\mathit{\eta }};q\left(t_n\right))\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{\eta }},\end{array}$$

while the X_i process considered in isolation using Eq. (4) results in an approximation to the average process rate of

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\underset{t_n}{\overset{t_n+\mathrm{\Delta }t}{\int }}{X}_i\left(q_{X_i}\right(\mathit{\eta }\left)\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\underset{t_n}{\overset{t_n+\mathrm{\Delta }t}{\int }}{X}_i\left(q_{X_i}\right[\mathit{\eta }-t_n;q_{X_i}^{\mathrm{input}}\left]\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}{X}_i\left(q_{X_i}\right[\stackrel{\mathrm{̃}}{\mathit{\eta }};q_{X_i}^{\mathrm{input}}\left]\right)\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{\eta }}.\end{array}$$

The discrepancy between these two average process rates has two sources:

1. The function $q_{X_i}$ differs from the function q, as the time evolution of $q_{X_i}$ is controlled by a single process (see Eq. 3), whereas the evolution of q is controlled by multiple processes (see Eq. 1).

2. The input state $q_{X_i}^{\mathrm{input}}$ used for integrating the equation of process X_i can differ from q(t_n).

In other words, the error in the estimated process rate X_i (or increment X_iΔt) can arise from (1) treating the process in isolation without considering the influence of other processes on the physical quantity, q, and (2) starting the single-process integration with an input that deviates from the solution of the multi-process ODE. These two types of error are referred to as isolation-induced error and input-induced error, respectively, in the remainder of the paper. To further elaborate on these two sources of splitting error, consider the following generic two-process ODE:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}q\left(t\right)}{\mathrm{d}t}}=A\left(q\right(t\left)\right)+B\left(q\right(t\left)\right),t>\mathrm{0},q\left(\mathrm{0}\right)=q^{\mathrm{IC}}.\end{array}$$

It has the following exact solution at t_n+1, written in terms of the exact solution at t_n:

$$q\left(t_{n+\mathrm{1}}\right)=q\left(t_n\right)+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right]\mathrm{d}\mathit{\eta }.$$

Below, the errors of two widely used coupling schemes are analyzed: parallel and sequential splitting. Both of these schemes are depicted in Fig. 1 using a flowchart description, a pseudo-code description, and an ODE description.

Figure 1 The parallel splitting (a1, b1, c1) and sequential splitting (a2, b2, c2) methods for solving the two-process ODE defined in Eq. (7). Panels (a1) and (a2), (b1) and (b2), and (c1) and (c2) depict the methods in three different ways.

2.2.1 Sources of error in parallel splitting

A parallel splitting scheme first lets each model component estimate the process rate of a single process and then sums up the corresponding increments to advance an input value of q to an output value of q. The scheme can be represented by a mapping $q^{\mathrm{output}}={\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{PS}}\left(q^{\mathrm{input}}\right)$, where

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{PS}}\left(q^{\mathrm{input}}\right)\equiv q^{\mathrm{input}}+\mathrm{\Delta }t\left(A^{*}+B^{*}\right),\\ {\displaystyle }& {\displaystyle }{A}^{*}\equiv \left(q_A\right[\mathrm{\Delta }t;q^{\mathrm{input}}]-q^{\mathrm{input}})/\mathrm{\Delta }t,\\ {\displaystyle }& {\displaystyle }{B}^{*}\equiv \left(q_B\right[\mathrm{\Delta }t;q^{\mathrm{input}}]-q^{\mathrm{input}})/\mathrm{\Delta }t.\end{array}$$

This gives

$$\begin{array}{ll}{\displaystyle }{\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{PS}}\left(q^{\mathrm{input}}\right)=q^{\mathrm{input}}& {\displaystyle }+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q_A\right[\mathit{\eta };q^{\mathrm{input}}\left]\right)\mathrm{d}\mathit{\eta }\\ \text{(8)}& {\displaystyle }& {\displaystyle }+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q^{\mathrm{input}}\left]\right)\mathrm{d}\mathit{\eta }.\end{array}$$

Recall that the local truncation error for a numerical method is, per the definition given in Eq. (2), the difference between the exact solution at t_n+1 and the solution at t_n+1 obtained from the numerical method applied to the exact solution at t_n. Thus, the local truncation error for parallel splitting is the difference between the exact two-process solution given by Eq. (5) and the result of Eq. (8) with q^input=q(t_n):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{PS}}\left(q\right(t_n\left)\right)-q\left(t_{n+\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }=\underset{{\mathrm{lte}}_A^{\mathrm{PS}}}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }}}\\ \text{(9)}& {\displaystyle }& {\displaystyle }+\underset{{\mathrm{lte}}_B^{\mathrm{PS}}}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }}}.\end{array}$$

Note that, as the time integration for A and B both start from q(t_n), there are only local truncation errors caused by treating A and B in isolation (i.e., no input-induced error). Also note that, as parallel splitting treats A and B the same way, the local truncation error expression shows symmetry between the two processes.

2.2.2 Sources of error in sequential splitting

The sequential splitting scheme discussed here handles different processes in a successive manner, letting each model compartment operate on an input value of q and return an updated value of q. Here, operator A is evaluated first, and the result is then used as input to operator B. The method can be represented by a mapping $q^{\mathrm{output}}={\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{SS}}\left(q^{\mathrm{input}}\right)$, where

$$\begin{array}{ll}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{SS}}\left(q^{\mathrm{input}}\right)\equiv q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q^{\mathrm{input}}\left]\right]{\displaystyle }& {\displaystyle }=q_A[\mathrm{\Delta }t;q^{\mathrm{input}}]+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q^{\mathrm{input}}]\left]\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }=q^{\mathrm{input}}+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q_A\right[\mathit{\eta };q^{\mathrm{input}}\left]\right)\mathrm{d}\mathit{\eta }\\ \text{(11)}& {\displaystyle }& {\displaystyle }+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q^{\mathrm{input}}]\left]\right)\mathrm{d}\mathit{\eta }.\end{array}$$

The local truncation error for sequential splitting is the difference between the exact two-process solution given by Eq. (5) and the result of Eq. (11) with q^input=q(t_n):

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& {\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{SS}}\left(q\right(t_n\left)\right)-q\left(t_{n+\mathrm{1}}\right)\\ & =\underset{{\mathrm{lte}}_A^{\mathrm{SS}}}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }}}\\ & +\underset{{\mathrm{lte}}_B^{\mathrm{SS}}}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q\right[\mathit{\eta },q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }}}.\end{array}\end{array}$$

Here, the ${\mathrm{lte}}_A^{\mathrm{SS}}$ term is the same as ${\mathrm{lte}}_A^{\mathrm{PS}}$ in parallel splitting (see Eq. 9) that includes only an isolation-induced error (q_A≠q), while the ${\mathrm{lte}}_B^{\mathrm{SS}}$ term includes not only an isolation-induced error (q_B≠q) but also an error resulting from the B process being integrated from an input that has already been updated by the A process, namely, q_A[Δt;q(t_n)].

2.3 The leading-order error terms

The analysis in the previous subsections provides a qualitative understanding of the sources of errors associated with different coupling methods. In order to obtain a more quantitative assessment of the error magnitudes and identify possible cancelation of different error terms, Taylor series expansion is used to derive the leading-order error terms of the local truncation error. The gist of the method is to expand the integrals in Sect. 2.2 about Δt=0, noting the following:

• $q_{X_i}\left(t\right)$ and q(t) are different functions whose time derivatives are given by Eqs. (1) and (3), respectively.

• The input state $q_{{}_{X_i}}^{\mathrm{input}}$ used for integrating the ODE of process X_i might deviate from the exact solution at t_n in a way that depends on Δt, and hence also need an expansion. For example, in the case of sequential splitting, $q_B^{\mathrm{input}}=q_A[\mathrm{\Delta }t;q(t_n\left)\right]$.

Appendix A explains in detail how the various integrals in Sect. 2.2 can be expanded to derive the leading-order error terms for the parallel and sequential splitting methods. Below, the key pieces of information obtained through these derivations are highlighted.

2.3.1 Leading-order errors in parallel splitting

The derivation detailed in Appendix A2 indicates that, when parallel splitting is used, the local truncation error attributable to the integration of process A (i.e., the part marked with ${\mathrm{lte}}_A^{\mathrm{PS}}$ in Eq. 9) is

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathrm{lte}}_A^{\mathrm{PS}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left(-{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}B\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

From the derivation in Appendix A2, it can be seen that the leading-order error (the (Δt)² term) results from the fact that the equation of $\mathrm{d}{q}_A/\mathrm{d}t$ lacks the B term that appears in the equation of $\mathrm{d}q/\mathrm{d}t$. In other words, the leading-order error is caused by integrating the A process in isolation without considering the influence of B. Similarly, the local truncation error attributable to the integration of process B (i.e., the part marked with ${\mathrm{lte}}_B^{\mathrm{PS}}$ in Eq. 9) is

$$\begin{array}{}\text{(14)}& {\displaystyle }{\mathrm{lte}}_B^{\mathrm{PS}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left(-{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}A\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

The leading-order error term is caused by applying the B process in isolation without considering the influence of A. The overall local truncation error for parallel splitting reads as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{PS}}\left(q\right(t_n\left)\right)-q\left(t_{n+\mathrm{1}}\right)={\mathrm{lte}}_A^{\mathrm{PS}}+{\mathrm{lte}}_B^{\mathrm{PS}}\\ \text{(15)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left(-{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}B-{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}A\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

As mentioned earlier, the expression has symmetry between the processes A and B, as the parallel splitting method treats the two processes in the same way.

2.3.2 Leading-order errors in sequential splitting

The derivation detailed in Appendix A3 shows that, when sequential splitting is used, the local truncation errors attributable to the integration of the respective A and B processes are as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }{\mathrm{lte}}_A^{\mathrm{SS}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left(-{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}B\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right),\text{(17)}& {\displaystyle }& {\displaystyle }{\mathrm{lte}}_B^{\mathrm{SS}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}(+{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}A){|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

Here, ${\mathrm{lte}}_A^{\mathrm{SS}}$ has the same expression as in parallel splitting (see Eq. 16 vs. Eq. 13); ${\mathrm{lte}}_B^{\mathrm{SS}}$ has the same form as in parallel splitting but a different sign, which results from the fact that the B process is integrated with an input state already updated by the A process, and the input-induced error overcompensates for the error caused by ignoring the influence of A when integrating the B process (see Appendix A3). The overall local truncation error for sequential splitting is as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{SS}}\left(q\right(t_n\left)\right)-q\left(t_{n+\mathrm{1}}\right)={\mathrm{lte}}_A^{\mathrm{SS}}+{\mathrm{lte}}_B^{\mathrm{SS}}\\ \text{(18)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left(-{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}B+{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}A\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

2.4 Framework summary and generalization

Leveraging the framework to derive and compare the splitting truncation errors of parallel and sequential splitting provides the following understanding that generalizes beyond the two-process problem. Note that a term containing $\left(\mathrm{d}B/\mathrm{d}q\right)A$ indicates an error caused by inaccurate accounts of the influence of process A on process B. This can be seen from the following Taylor expansion:

$$\begin{array}{ll}{\displaystyle }\mathrm{\Delta }t\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}A\right){|}_{q=q\left(t_n\right)}=& {\displaystyle }B\left(q\right(t_n)+\mathrm{\Delta }tA(q\left(t_n\right)\left)\right)-B\left(q\right(t_n\left)\right)\\ {\displaystyle }& {\displaystyle }+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{2}}\right).\end{array}$$

A negative sign in front of $\left(\mathrm{d}B/\mathrm{d}q\right)A$ suggests a lack or underestimation of the influence of A on B, whereas a positive sign suggests an overestimation of that influence. Isolation-induced error in the form of integrating B in isolation will lead to an underestimation of the influence of A on B. If there are additional processes, integrating B in isolation will lead to an underestimation of the influence of all other processes on B. Input-induced error in the form of using an input state updated by a full time step worth of A will lead to an overestimation of the influence of A on B. If there are additional processes, using an input state updated by a full time step worth of A and other processes will lead to an overestimation of A and those other processes on B. With this understanding, the framework is easily generalized to coupling methods that utilize different combinations of input states across numerous processes. The following section will demonstrate such a generalization to a three-process problem.

It is also worthwhile to note that the framework presented here can be revised to evaluate the overall truncation error in a temporally discretized system. Namely, one would replace $q_{X_i}[t-t_n;q_{X_i}^{\mathrm{input}}]$ in Appendix A1 with the time integration scheme used for process X_i. For example, if forward Euler is used for integration of process X_i, one would use

$${\displaystyle }{q}_{X_i}[t-t_n;q_{X_i}^{\mathrm{input}}]\equiv q_{X_i}^{\mathrm{input}}+(t-t_n)X_i\left(q_{X_i}^{\mathrm{input}}\right).$$

Such a revised framework would lead to results equivalent to those presented in studies such as Ubbiali et al. (2021). That said, it is also worthwhile to note that such overall truncation error results can be more difficult to interpret than the results of the framework presented here. Take, for example, how the overall truncation error for parallel splitting with forward Euler time integration of all processes is identical to the overall truncation for an unsplit forward Euler approach, which may seem to suggest that there is no splitting truncation error for the parallel splitting method in general. However, the overall truncation error expressions for parallel splitting and the unsplit approach will not be identical if backward Euler time integration is instead used, revealing that there is indeed splitting truncation error for the parallel splitting method.

3 The semi-discrete analysis framework applied to a three-process problem inspired by EAMv1

The semi-discrete error analysis framework presented in Sect. 2 is now used to analyze the dust life cycle problem in EAMv1. The dust mass budget analyses carried out using EAMv1's simulation output and presented in Sect. 3.1 of Part 1 (Wan et al., 2024) have revealed that, at the global scale, the strongest sources and sinks of dust aerosols are (i) surface emissions, (ii) dry removal, and (iii) turbulent mixing and aerosol activation–resuspension. As such, we focus on these sources and sinks, ignore the many other aerosol-related processes in EAM, and consider a canonical three-process problem:

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\mathrm{d}q\left(t\right)}{\mathrm{d}t}}=A\left(q\right(t\left)\right)+B\left(q\right(t\left)\right)+C\left(q\right(t\left)\right),t>\mathrm{0},q\left(\mathrm{0}\right)=q^{\mathrm{IC}},\end{array}$$

where q is a dust mass mixing ratio, A represents the emissions, B represents dry removal, and C corresponds to turbulent mixing. As in Sect. 2, denote a discrete time step Δt and discrete time points $t_{n+\mathrm{1}}=t_n+\mathrm{\Delta }t$. Denote the numerical solution at time t_n+1 as qⁿ⁺¹ and the exact solution it approximates as q(t_n+1). Figure 2 describes two schemes for obtaining qⁿ⁺¹ from qⁿ, corresponding to the original and revised process coupling schemes in EAMv1 discussed in Wan et al. (2024).

Figure 2 Three different descriptions of two process coupling schemes for the three-process problem defined in Sect. 3. The scheme depicted in panels (a1), (b1), and (c1) corresponds to the original scheme used in EAMv1 for the coupling of aerosol emissions, dry removal, and the parameterization of turbulent transport and aerosol activation–resuspension. The scheme depicted in the panels (a2), (b2), and (c2) corresponds to the revised scheme proposed and evaluated in Part 1 (Wan et al., 2024). We note that these descriptions are simplified versions of the coupling implemented in EAMv1. Here, we focus only on the three strongest sources and sinks of the global mean dust budget presented in Sect. 3 of Wan et al. (2024), while the many other processes in EAMv1 (see Fig. 1 in Wan et al., 2024) are omitted.

Consider three single-process ODEs in the form of Eq. (3) where $X_i=A,B,$ or C, namely,

$$\begin{array}{l}{\displaystyle \frac{\mathrm{d}{q}_A\left(t\right)}{\mathrm{d}t}}=A\left(q\right(t\left)\right),t>t_n,q_A\left(t_n\right)=q_A^{\mathrm{input}},\\ {\displaystyle \frac{\mathrm{d}{q}_B\left(t\right)}{\mathrm{d}t}}=B\left(q\right(t\left)\right),t>t_n,q_B\left(t_n\right)=q_B^{\mathrm{input}},\\ {\displaystyle \frac{\mathrm{d}{q}_C\left(t\right)}{\mathrm{d}t}}=C\left(q\right(t\left)\right),t>t_n,q_C\left(t_n\right)=q_C^{\mathrm{input}}.\end{array}$$

The exact solutions are denoted using the notation defined in Eq. (4), namely,

$$\begin{array}{ll}{\displaystyle }{q}_A[t-t_n;q_A^{\mathrm{input}}]& {\displaystyle }\equiv q_A^{\mathrm{input}}+\underset{t_n}{\overset{t}{\int }}A\left(q_A\right(\mathit{\eta }\left)\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }=q_A^{\mathrm{input}}+\underset{\mathrm{0}}{\overset{t-t_n}{\int }}A\left(q_A\right[\stackrel{\mathrm{̃}}{\mathit{\eta }};q_A^{\mathrm{input}}\left]\right)\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{\eta }},\end{array}$$

$$\begin{array}{ll}{\displaystyle }{q}_B[t-t_n;q_B^{\mathrm{input}}]& {\displaystyle }\equiv q_B^{\mathrm{input}}+\underset{t_n}{\overset{t}{\int }}B\left(q_B\right(\mathit{\eta }\left)\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }=q_B^{\mathrm{input}}+\underset{\mathrm{0}}{\overset{t-t_n}{\int }}B\left(q_B\right[\stackrel{\mathrm{̃}}{\mathit{\eta }};q_B^{\mathrm{input}}\left]\right)\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{\eta }},\end{array}$$

$$\begin{array}{ll}{\displaystyle }{q}_C[t-t_n;q_C^{\mathrm{input}}]& {\displaystyle }\equiv q_C^{\mathrm{input}}+\underset{t_n}{\overset{t}{\int }}C\left(q_C\right(\mathit{\eta }\left)\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }=q_C^{\mathrm{input}}+\underset{\mathrm{0}}{\overset{t-t_n}{\int }}C\left(q_C\right[\stackrel{\mathrm{̃}}{\mathit{\eta }};q_C^{\mathrm{input}}\left]\right)\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{\eta }}.\end{array}$$

3.1 Process coupling schemes

The original EAMv1 uses sequential splitting for all three processes (see Fig. 2a1, b1, c1), which can be represented by the mapping

$$\begin{array}{}\text{(20)}& {\displaystyle }{q}^{n+\mathrm{1}}={\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{Ori}}\left(q^n\right)\equiv q_C[\mathrm{\Delta }t;q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q^n]\left]\right]\end{array}$$

or, equivalently,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{q}^{n+\mathrm{1}}=q_C[\mathrm{\Delta }t;q^n+(q_A[\mathrm{\Delta }t;q^n]-q^n)\\ \text{(21)}& {\displaystyle }& {\displaystyle }+\left(q_B\right[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q^n]]-q_A[\mathrm{\Delta }t;q^n\left]\right)].\end{array}$$

Defining

$$\begin{array}{}\text{(22)}& {\displaystyle }{A}^{*}\equiv \left(q_A\right[\mathrm{\Delta }t;q^n]-q^n)/\mathrm{\Delta }t,\text{(23)}& {\displaystyle }{B}^{*}\equiv \left(q_B\right[\mathrm{\Delta }t;q^n]-q^n)/\mathrm{\Delta }t,\end{array}$$

the revised coupling scheme depicted in Fig. 2a2, b2, and c2 can be represented by the mapping

$${\displaystyle }{q}^{n+\mathrm{1}}={\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{Rev}}\left(q^n\right)\equiv q_C[\mathrm{\Delta }t;q^n+\mathrm{\Delta }t(A^{*}+B^{*}\left)\right]$$

or, equivalently,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{q}^{n+\mathrm{1}}=q_C[\mathrm{\Delta }t;q^n+(q_A[\mathrm{\Delta }t;q^n]-q^n)\\ \text{(24)}& {\displaystyle }& {\displaystyle }+\left(q_B\right[\mathrm{\Delta }t;q^n]-q^n)].\end{array}$$

3.2 Leading-order error terms

The original and revised coupling schemes described in Eqs. (20) and (24) can be viewed as different combinations of the two-process sequential and parallel splitting schemes discussed in Sect. 2. Based on the discussions in that section, and keeping in mind that the focus here is the local truncation error, one can make the following reasoning about the original coupling scheme in EAMv1:

• For process A, as the solution procedure starts from the exact solution at t_n and integrates the A term in isolation, one expects to get order (Δt)² errors caused by performing time integration without considering the impacts of B and C on the A process. The coefficients in front of $\left(\mathrm{d}A/\mathrm{d}q\right)B$ and $\left(\mathrm{d}A/\mathrm{d}q\right)C$ are expected to be $-(\mathrm{\Delta }t{)}^{\mathrm{2}}/\mathrm{2}$. In other words, the splitting truncation error associated with the A process is expected to be

  $$\begin{array}{}\text{(25)}& {\mathrm{lte}}_A^{\mathrm{Ori}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}(-B-C)\right]{|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

• For process B, as the solution procedure starts from a mixing ratio updated by A and ignores the A and C terms on the RHS of the original ODE, one expects there to be two error terms caused by integrating B in isolation and an error associated with the input state. The input-induced error is expected to overcompensate for the error caused by ignoring the impact of A on the B process. Hence, the splitting truncation error associated with the B process is expected to have the form

  $$\begin{array}{}\text{(26)}& {\mathrm{lte}}_B^{\mathrm{Ori}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}(+A-C)\right]{|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

• For process C, the solution procedure starts from a mixing ratio updated by both A and B, and the C term is integrated in isolation. Therefore, one expects to have two input-induced errors overcompensating for two isolation-induced errors, giving a splitting truncation error in the form of

  $$\begin{array}{}\text{(27)}& {\mathrm{lte}}_C^{\mathrm{Ori}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}(+A+B)\right]{|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

The overall local truncation error in the original coupling is expected to be

$$\begin{array}{}\text{(28)}& \begin{array}{rl}& {\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{Ori}}\left(q\right(t_n\left)\right)-q\left(t_{n+\mathrm{1}}\right)={\mathrm{lte}}_A^{\mathrm{Ori}}+{\mathrm{lte}}_B^{\mathrm{Ori}}+{\mathrm{lte}}_C^{\mathrm{Ori}}\\ & ={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}(-B-C)+{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}(A-C)+{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}(A+B)\right]{|}_{q=q\left(t_n\right)}\\ & +\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}\end{array}$$

The revised coupling scheme differs from the original scheme only in the input state for the integration of the B process; see Eq. (24) vs. Eq. (21). Therefore, one expects that the splitting truncation errors associated with the other two processes, A and C, are the same as in the original scheme and that the splitting truncation error associated with the B process has a minus sign instead of a plus sign for the $\left(\mathrm{d}B/\mathrm{d}q\right)A$ term (i.e., no input-induced error, only the isolation-induced error). In other words,

$$\begin{array}{}\text{(29)}& {\displaystyle }{\mathrm{lte}}_A^{\mathrm{Rev}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}(-B-C)\right]{|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right),\end{array}$$

$$\begin{array}{}\text{(30)}& {\displaystyle }{\mathrm{lte}}_B^{\mathrm{Rev}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}(-A-C)\right]{|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right),\end{array}$$

$$\begin{array}{}\text{(31)}& {\displaystyle }{\mathrm{lte}}_C^{\mathrm{Rev}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}(+A+B)\right]{|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right),\end{array}$$

and the overall local truncation error is expected to be

$$\begin{array}{}\text{(32)}& \begin{array}{rl}& {\mathcal{F}}_{\mathrm{\Delta }t}^{\mathrm{Rev}}\left(q\right(t_n\left)\right)-q\left(t_{n+\mathrm{1}}\right)={\mathrm{lte}}_A^{\mathrm{Rev}}+{\mathrm{lte}}_B^{\mathrm{Rev}}+{\mathrm{lte}}_C^{\mathrm{Rev}}\\ & ={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}(-B-C)+{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}(-A-C)+{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}(A+B)\right]{|}_{q=q\left(t_n\right)}\\ & +\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}\end{array}$$

All of the error expressions in Eqs. (25)–(32) are confirmed by the step-by-step derivations presented in Appendix B. This agreement demonstrates how the two-process splitting truncation error results can be leveraged as building blocks to derive splitting truncation errors for multi-process problems using logical reasoning instead of Taylor series expansions. In other words, the derivation of splitting truncation error for multi-process problems does not always have to be done in the step-by-step manner in Appendix B. The logical reasoning approach not only can facilitate rapid development of alternative coupling schemes but also makes the framework more accessible to model developers on the physics side who might find the lengthy calculus derivations too tedious or daunting.

3.3 Characteristics of the leading-order error terms inferred from EAMv1 results

Recall that the three RHS terms in the three-process ODE discussed above are meant to represent the respective surface emissions (A), dry removal (B), and turbulent mixing (C) of dust aerosols in EAMv1. The parameterization descriptions and EAMv1 simulations presented in Sects. 2 and 3 of Part 1 (Wan et al., 2024) can be used to infer several features of the leading-order error terms listed above in Sect. (3.2). The dust budget analyses shown in Sect. 3 of Wan et al. (2024) indicate that the dominant sources and sinks are found in the lowest model layer in the dust source regions, where A (emission) is a source, B (dry removal) is typically a sink, and C (turbulent mixing) is typically a sink (i.e., A>0, B<0, and C<0). Given the same air density and deposition velocity, the downward dry removal flux at the Earth's surface is proportional to the mean dust mixing ratio of the layer (see Eq. 1 in Wan et al., 2024). This means that one can expect

$$\begin{array}{}\text{(33)}& {\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}<\mathrm{0}\end{array}$$

to be true in typical cases. Equation (33) is confirmed by the scatterplot in Fig. 3, where the dry removal rate B is plotted against dust mixing ratio q using 90 d of 6-hourly output in dust source regions simulated with the original EAMv1.

Figure 3 Dust aerosol dry removal rate (y axis) plotted against dust aerosol mixing ratio (x axis) in the lowest model layer in dust source regions simulated by the original EAMv1 using a vertical grid with 72 layers. The data used in the figure included 90 d of 6-hourly instantaneous output.

It then follows that the local truncation error associated with the B process in EAMv1's original process coupling (see Eq. 26) can be written as

$$\begin{array}{}\text{(34)}& {\mathrm{lte}}_B^{\mathrm{Ori}}=-{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left[\left|{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}\right|\left(\right|A|+|C\left|\right)\right]{|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right),\end{array}$$

and the local truncation error of the B process in the revised coupling scheme (see Eq. 30) can be written as

$$\begin{array}{}\text{(35)}& {\mathrm{lte}}_B^{\mathrm{Rev}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left[\left|{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}\right|\left(\right|A|-|C\left|\right)\right]{|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

Because $\left|\right|A|-|C\left|\right|\le \left|A\right|+\left|C\right|$, with equality only when A or C is zero, it is expected that the magnitude of the leading-order local truncation error associated with process B is smaller in the revised coupling than in the original scheme, i.e.,

$$\begin{array}{}\text{(36)}& {\displaystyle }\left|{\mathrm{lte}}_B^{\mathrm{Rev}}\right|\lesssim \left|{\mathrm{lte}}_B^{\mathrm{Ori}}\right|.\end{array}$$

This result, combined with the fact that the local truncation errors associated with the other two processes (A and C) have the same expressions in the original and revised coupling schemes, provides a justification for adopting the revised scheme in EAMv1.

To see that Eq. (36) holds in practice, it is useful to first note that the main leading-order difference between ${\mathrm{lte}}_B^{\mathrm{Ori}}$ in Eq. (26) and ${\mathrm{lte}}_B^{\mathrm{Rev}}$ in Eq. (30) is the term (A−C) vs. $(-A-C)$, respectively, evaluated at q=q(t_n). While the values of A(q(t_n)) and C(q(t_n)) are not feasible to obtain in practice, as q(t_n) is the exact solution, the approximate A and C values calculated and used in EAMv1 simulations are relatively straightforward to obtain using the online diagnostic tool of Wan et al. (2022), as was done in Part 1 (Wan et al., 2024). Figure 4 shows the annual mean values of the computed (A−C) and $(-A-C)$ in dust source regions in North Africa, using both the original and revised coupling methods.

Figure 4 Comparison of key terms in the splitting truncation error associated with dry removal (process B) using 10-year mean interstitial dust mass mixing ratio process rates (unit: kg kg⁻¹ s⁻¹) caused by emissions (process A) and dry removal (process B) in the lowest model layer in EAMv1 simulations using the original coupling scheme (a, b) and the revised scheme (c, d).

When the original coupling method is used, the magnitude of (A−C) ranges from being slightly to substantially larger than the magnitude of $(-A-C)$. When the revised coupling is used, the magnitude of (A−C) is substantially larger everywhere compared with that of $(-A-C)$. Both results support that the leading-order term in ${\mathrm{lte}}_B^{\mathrm{Rev}}$ is smaller in magnitude than that of the leading-order term in ${\mathrm{lte}}_B^{\mathrm{Ori}}$, i.e., Eq. (36).

It is also worth noting that the leading-order term in ${\mathrm{lte}}_B^{\mathrm{Ori}}$ is negative (by Eq. 34), and, as B itself is negative, the negative leading-order term in the truncation error indicates an overestimation of the B process at each time step of the original coupling method. The study by Feng et al. (2022) has pointed out that dust dry removal in EAMv1 is generally overestimated in dust source regions; thus, the reduction in |lte_B| shown in Eq. (36) is consistent with the significantly weaker dry removal seen in Part 1 (Wan et al., 2024) when the revised method is used instead of the original method in EAMv1. While the local truncation error caused by process splitting is not the only source of error in global simulations (other error sources include propagated splitting error in Eq. 2, temporal and spatial discretization errors in dry removal and other aerosol processes, model formulation error, and parameter uncertainty, among others), the theoretical analysis here and the global simulations in Feng et al. (2022) and Wan et al. (2024) suggest that ${\mathrm{lte}}_B^{\mathrm{Ori}}$ is likely an important contributor to the overly strong dust dry removal in the original EAMv1.

4 Summary and conclusions

A semi-discrete error analysis framework was introduced for assessing splitting error of process coupling methods. By assuming that the time integration of each individual process is exact, the framework identified two general sources of splitting error. The first is denoted isolation-induced error and is from the integration, exact or otherwise, of a process without the influence of other processes (i.e., in isolation). The second is denoted input-induced error and is from starting the single-process integration from an input that deviates from the full solution of the multi-process equation. The corresponding splitting truncation error terms from those two sources were derived for a generic two-process problem for two common coupling methods: parallel splitting and sequential splitting. The parallel splitting method results in isolation-induced error from both processes. The combination of isolation-induced errors leads to an underestimation of the influence of the one process on the other. The sequential splitting method results in isolation-induced error from the first process that is integrated in a time step and in both types of error for the second process. The combination of the two types of errors leads to an underestimation of the influence of the second process on the first, as in parallel splitting, but the input-induced error from the second process is shown to overcompensate for the isolation-induced error from the first process. Thus, sequential splitting results in an overestimate of the influence of the first process on the second process.

A three-process problem and two coupling schemes were analyzed that corresponded to the coupling of dust emissions, dry removal, and turbulence mixing in the original EAMv1 and the revised coupling scheme proposed in Part 1 (Wan et al., 2024). The semi-discrete analysis revealed that the original and revised coupling schemes have the same forms of leading-order error for the emissions and turbulent mixing, while the magnitude of the local truncation error in dry removal is smaller in the revised scheme. Assuming it is useful, especially in the long run, to address different error sources in EAM separately, this result provides a justification to adopt the revised coupling in EAMv1. The analysis also revealed that the local truncation error of the dry removal process in the original EAMv1 corresponds to an overestimation of the dry removal rate. This result, combined with the EAMv1 results presented in Feng et al. (2022) and Wan et al. (2024), suggests that the sequential splitting of emissions, dry removal, and turbulent mixing is likely an important contributor to the overestimated dry removal in dust source regions in the original EAMv1.

While the error analysis framework is presented in the context of discussions on the dust life cycle in EAMv1, using the framework as done in this work is much more general. For one, such a framework can be used to analyze coupling methods beyond the two schemes discussed here and in Part 1 (Wan et al., 2024) as well as numerical coupling problems involving more than three processes. Additionally, because many applications rely on low-order coupling methods such as those discussed here, this paper shows how such a framework can be used to inform coupling approach choices in areas beyond atmospheric climate, including hydrology, fusion, reactive flow modeling, and many others. As such, the authors plan to use the framework to help further reduce splitting errors in EAM and other applications.

Appendix A: Derivation of the leading-order error terms in a generic two-process ODE

This section details the step-by-step derivation of the leading-order local truncation error terms caused by applying the parallel splitting and sequential splitting methods to solving the two-process problem defined in Eq. (7). The starting point is the local truncation error terms ${\mathrm{lte}}_A^{\text{SS}}$ and ${\mathrm{lte}}_B^{\text{SS}}$ in Eq. (12) and ${\mathrm{lte}}_A^{\text{PS}}$ and ${\mathrm{lte}}_B^{\text{PS}}$ in Eq. (9).

A1 Taylor expansion of an integral

The Taylor expansion of an integral is a key element in deriving the leading-order terms. As such, it is formalized herein. Recall from Eq. (5) that q[η;ϕ] is the exact solution of the multi-process problem Eq. (1) evolved from input state ϕ for time η. For a function f(δ) defined as

$${\displaystyle }f\left(\mathit{\delta }\right)=\underset{\mathrm{0}}{\overset{\mathit{\delta }}{\int }}F\left(q\right[\mathit{\eta };\mathit{\varphi }\left]\right)\mathrm{d}\mathit{\eta },$$

the first and second derivatives are

$$\begin{array}{ll}{\displaystyle }{f}^{\prime }\left(\mathit{\delta }\right)& {\displaystyle }=F\left(q\right[\mathit{\delta };\mathit{\varphi }\left]\right)\text{and}\\ {\displaystyle }{f}^{\prime \prime }\left(\mathit{\delta }\right)& {\displaystyle }={\displaystyle \frac{\mathrm{d}F}{\mathrm{d}q}}\left(q\right[\mathit{\delta };\mathit{\varphi }\left]\right){\displaystyle \frac{\mathrm{d}q}{\mathrm{d}t}}[\mathit{\delta };\mathit{\varphi }],\end{array}$$

respectively, and thus

$$\begin{array}{ll}{\displaystyle }f\left(\mathrm{0}\right)& {\displaystyle }=\mathrm{0},\\ {\displaystyle }{f}^{\prime }\left(\mathrm{0}\right)& {\displaystyle }=F\left(q\right[\mathrm{0};\mathit{\varphi }\left]\right)=F\left(\mathit{\varphi }\right),\\ {\displaystyle }{f}^{\prime \prime }\left(\mathrm{0}\right)& {\displaystyle }={\displaystyle \frac{\mathrm{d}F}{\mathrm{d}q}}\left(q\right[\mathrm{0},\mathit{\varphi }\left]\right){\displaystyle \frac{\mathrm{d}q}{\mathrm{d}t}}[\mathrm{0};\mathit{\varphi }]={\displaystyle \frac{\mathrm{d}F}{\mathrm{d}q}}\left(\mathit{\varphi }\right){\displaystyle \frac{\mathrm{d}q}{\mathrm{d}t}}[\mathrm{0};\mathit{\varphi }].\end{array}$$

Note that Eq. (6) was used to simplify the expression. The Taylor expansion of f(Δt) about Δt=0 is now given as

$$\begin{array}{ll}{\displaystyle }f\left(\mathrm{\Delta }t\right)=& {\displaystyle }f\left(\mathrm{0}\right)+\mathrm{\Delta }t{f}^{\prime }\left(\mathrm{0}\right)+{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}{f}^{\prime \prime }\left(\mathrm{0}\right)+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right)\\ {\displaystyle }=& {\displaystyle }\mathrm{0}+\mathrm{\Delta }tF\left(\mathit{\varphi }\right)+{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{d}F}{\mathrm{d}q}}\left(\mathit{\varphi }\right){\displaystyle \frac{\mathrm{d}q}{\mathrm{d}t}}[\mathrm{0};\mathit{\varphi }]\\ \text{(A1)}& {\displaystyle }& {\displaystyle }+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

Note that Eq. (A1) can be used to show both that

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& A\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tA\left(q\right(t_n\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}(A+B)\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}})\end{array}\end{array}$$

and

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& B\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tB\left(q\right(t_n\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}(A+B)\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}\end{array}$$

A2 Parallel splitting

Recall the local truncation error term for parallel splitting in Eq. (9):

$$\begin{array}{}\text{(A4)}& {\displaystyle }{\mathrm{lte}}_A^{\text{PS}}=\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }.\end{array}$$

The second integral is expanded using Eq. (A2). For the first integral, use Eq. (A1), with F=A, q=q_A, and ϕ=q(t_n) so that $F\left(q\right[\mathit{\eta };\mathit{\varphi }\left]\right)=A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)$, to find

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tA\left(q\right(t_n\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}\left(q\right(t_n\left)\right){\displaystyle \frac{\mathrm{d}{q}_A}{\mathrm{d}t}}[\mathrm{0};q(t_n\left)\right]+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),\end{array}$$

which can be simplified using Eq. (3) to get

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tA\left(q\right(t_n\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}A\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}$$

The expansions of the integrals in Eq. (A4) are now combined to find

$${\displaystyle }{\mathrm{lte}}_A^{\text{PS}}=-{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}B\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),$$

as shown in Eq. (13). Recall the other local truncation error term from Eq. (9):

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\mathrm{lte}}_B^{\text{PS}}=\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }.\end{array}$$

The second integral is expanded using Eq. (A3). For the first integral, use Eq. (A1), with F=B, q=q_B, and ϕ=q(t_n) so that $F\left(q\right[\mathit{\eta };\mathit{\varphi }\left]\right)=B\left(q_B\right[\mathit{\eta };q\left(t_n\right)\left]\right)$, to find

$$\begin{array}{}\text{(A6)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{\int }_{\mathrm{0}}^{\mathrm{\Delta }t}& {\displaystyle }B\left(q_B\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tB\left(q\right(t_n\left)\right)\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}\left(q\right(t_n\left)\right){\displaystyle \frac{\mathrm{d}{q}_B}{\mathrm{d}t}}[\mathrm{0};q(t_n\left)\right]\\ {\displaystyle }& {\displaystyle }+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right),\end{array}\end{array}$$

which can be simplified using Eq. (3) to get

$$\begin{array}{ll}{\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& {\displaystyle }B\left(q_B\right[\mathit{\eta };q^n\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tB\left(q\right(t_n\left)\right)\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}B\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

The expansions of the integrals in Eq. (A5) are now combined to find

$${\displaystyle }{\mathrm{lte}}_B^{\text{PS}}=-{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}A\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),$$

as shown in Eq. (14).

A3 Sequential splitting

Recall the local truncation error term for sequential splitting in Eq. (12):

$${\displaystyle }{\mathrm{lte}}_A^{\text{SS}}=\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }.$$

Note that ${\mathrm{lte}}_A^{\text{SS}}$ is equivalent to ${\mathrm{lte}}_A^{\text{PS}}$, which has already been derived in Sect. A2 and is equivalent to Eq. (16). Recall the other local truncation error term from Eq. (12):

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}{\mathrm{lte}}_B^{\text{SS}}& =\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\mathrm{d}\mathit{\eta }\\ & -\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }.\end{array}\end{array}$$

The second integral is expanded using Eq. (A3). For the first integral, use Eq. (A1), with F=B, q=q_B, and $\mathit{\varphi }=q_A[\mathrm{\Delta }t;q(t_n\left)\right]$ so that $F\left(q\right[\mathit{\eta };\mathit{\varphi }\left]\right)=B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right])$, to find

$$\begin{array}{ll}{\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& {\displaystyle }B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tB\left(q_A\right[\mathrm{\Delta }t;q\left(t_n\right)\left]\right])\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}\left(q_A\right[\mathrm{\Delta }t;q\left(t_n\right)\left]\right){\displaystyle \frac{\mathrm{d}{q}_B}{\mathrm{d}t}}[\mathrm{0},q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]\\ {\displaystyle }& {\displaystyle }+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right),\end{array}$$

which can be simplified using Eq. (3) to get

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tB\left(q_A\right[\mathrm{\Delta }t;q\left(t_n\right)\left]\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}B\right){|}_{q=q_A[\mathrm{\Delta }t;q(t_n\left)\right]}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}\end{array}$$

To continue the expansion, use

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{q}_A[\mathrm{\Delta }t;q(t_n\left)\right]& =q_A[\mathrm{0};q(t_n\left)\right]+\mathrm{\Delta }t{\displaystyle \frac{\mathrm{d}{q}_A}{\mathrm{d}t}}[\mathrm{0};q(t_n\left)\right]+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{2}}\right)\\ & =q\left(t_n\right)+\mathrm{\Delta }tA\left(q\right(t_n\left)\right)+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{2}}\right)\end{array}\end{array}$$

to get

$$\begin{array}{ll}{\displaystyle }\mathrm{\Delta }tB& {\displaystyle }\left(q_A\right[\mathrm{\Delta }t;q\left(t_n\right)\left]\right)\\ {\displaystyle }& {\displaystyle }=\mathrm{\Delta }tB\left(q\right(t_n\left)\right)+\left(\mathrm{\Delta }t{)}^{\mathrm{2}}\right({\displaystyle \frac{dB}{\mathrm{d}q}}A){|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right)\end{array}$$

and

$$\begin{array}{ll}{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}& {\displaystyle }\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}B\right){|}_{q=q_A[\mathrm{\Delta }t;q^n]}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}B\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

This allows for the simplification of Eq. (A8) to

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\mathrm{d}\mathit{\eta }\\ & =\mathrm{\Delta }tB\left(q\right(t_n\left)\right)+\left(\mathrm{\Delta }t{)}^{\mathrm{2}}\right({\displaystyle \frac{dB}{\mathrm{d}q}}A){|}_{q=q\left(t_n\right)}+{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}B){|}_{q=q\left(t_n\right)}\\ & +\mathcal{O}\left((\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}\end{array}$$

The expansions of the integrals in Eq. (A7) are now combined to find

$${\displaystyle }{\mathrm{lte}}_B^{\text{SS}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}A\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),$$

as shown in Eq. (17). It is worth noting that compared with the expansion of Eq. (A6), the expansion of Eq. (A8) has an extra term $\left(\mathrm{\Delta }t{)}^{\mathrm{2}}\right(\frac{\mathrm{d}B}{\mathrm{d}q}A){|}_{q=q\left(t_n\right)}$ that results from the sequential splitting method using q_A[Δt;q(t_n)] (the value of q already updated by process A) as the input when integrating the $\mathrm{d}{q}_B/\mathrm{d}t$ equation. This leads to the sign difference between ${\mathrm{lte}}_B^{\text{SS}}$ and ${\mathrm{lte}}_B^{\text{PS}}$ that can be traced to the fact that the input used when integrating process B, i.e., q_A[Δt;q(t_n)], results in a leading-order error in the solution that overcompensates for the leading-order term caused by integrating the equation of $\mathrm{d}{q}_B/\mathrm{d}t$ without an A term on the RHS.

Appendix B: Derivation of the leading-order error terms in a three-process ODE

This section details the derivation of the leading-order local truncation error terms caused by the original splitting and revised splitting methods to solving the three-process problem defined in Eq. (19). As it will be useful herein, start by using Eq. (A1) to show that

$$\begin{array}{}\text{(B1)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& A\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tA\left(q\right(t_n\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}(A+B+C)\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),\end{array}\text{(B2)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& B\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tB\left(q\right(t_n\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}(A+B+C)\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),\end{array}\end{array}$$

and

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& C\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tC\left(q\right(t_n\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}(A+B+C)\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}\end{array}$$

B1 Revised splitting

The mapping in the revised splitting described in Eq. (24) can be written as

$$\begin{array}{rl}{\mathcal{F}}_{\mathrm{\Delta }t}^{\text{Rev}}\left(q\right(t_n\left)\right)& =q_C[\mathrm{\Delta }t;q(t_n)+\mathrm{\Delta }t(A^{*}+B^{*}\left)\right]\\ & =q\left(t_n\right)+\mathrm{\Delta }t(A^{*}+B^{*})\\ & +\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}C\left(q_C\right[\mathit{\eta };q\left(t_n\right)+\mathrm{\Delta }t(A^{*}+B^{*})\left]\right)\mathrm{d}\mathit{\eta }\\ & =q\left(t_n\right)+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }\\ & +\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }\\ & +\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}C\left(q_C\right[\mathit{\eta };q\left(t_n\right)+\mathrm{\Delta }t(A^{*}+B^{*})\left]\right)\mathrm{d}\mathit{\eta }.\end{array}$$

Thus, the local truncation error is expressed as

$$\begin{array}{rl}{\mathcal{F}}_{\mathrm{\Delta }t}^{\text{Rev}}\left(q\right(t_n\left)\right)& -q\left(t_{n+\mathrm{1}}\right)=\underset{{\mathrm{lte}}_A^{\text{Rev}}}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }}}\\ & +\underset{{\mathrm{lte}}_B^{\text{Rev}}}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }}}\\ & +\underset{{\mathrm{lte}}_C^{\text{Rev}}}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}C\left(q_C\right[\mathit{\eta };q\left(t_n\right)+\mathrm{\Delta }t(A^{*}+B^{*})\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}C\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }}}.\end{array}$$

The leading-order terms in ${\mathrm{lte}}_A^{\text{Rev}}$, ${\mathrm{lte}}_B^{\text{Rev}}$, and ${\mathrm{lte}}_C^{\text{Rev}}$ will now be derived in a manner similar to their two-process counterparts in Sect. A. For ${\mathrm{lte}}_A^{\text{Rev}}$, the second integral is expanded using Eq. (B1). For the first integral in ${\mathrm{lte}}_A^{\text{Rev}}$, use Eq. (A1), with F=A, q=q_A, and ϕ=q(t_n) so that $F\left(q\right[\mathit{\eta };\mathit{\varphi }\left]\right)=A\left(q_A\right[\mathit{\eta };\mathit{\varphi }\left]\right)$, to find

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tA\left(q\right(t_n\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}\left(q\right(t_n\left)\right){\displaystyle \frac{\mathrm{d}{q}_A}{\mathrm{d}t}}[\mathrm{0};q(t_n\left)\right]+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),\end{array}$$

which can be simplified using Eq. (6) to get

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tA\left(q\right(t_n\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}A\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}\end{array}$$

The expansions of the integrals in ${\mathrm{lte}}_A^{\text{Rev}}$ are now combined to find

$${\displaystyle }{\mathrm{lte}}_A^{\text{Rev}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left(-{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}q}}(B+C)\right)+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),$$

which is equivalent to the expression in Eq. (29). For ${\mathrm{lte}}_B^{\text{Ref}}$, the second integral is expanded using Eq. (B2). For the first integral in ${\mathrm{lte}}_B^{\text{Rev}}$, use Eq. (A1), with F=B, q=q_B, and ϕ=q(t_n) so that $F\left(q\right[\mathit{\eta };\mathit{\varphi }\left]\right)=B\left(q_B\right[\mathit{\eta };\mathit{\varphi }\left]\right)$, to find

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& B\left(q_B\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tB\left(q\right(t_n\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}\left(q\right(t_n\left)\right){\displaystyle \frac{\mathrm{d}{q}_B}{\mathrm{d}t}}[\mathrm{0};q(t_n\left)\right]+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),\end{array}$$

which can be simplified using Eq. (6) to get

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& B\left(q_B\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tB\left(q\right(t_n\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}B\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}\end{array}$$

The expansions of the integrals in ${\mathrm{lte}}_B^{\text{Rev}}$ are now combined to find

$${\displaystyle }{\mathrm{lte}}_B^{\text{Rev}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left(-{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}(A+C)\right)+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),$$

which is equivalent to the expression in Eq. (30). For ${\mathrm{lte}}_C^{\text{Rev}}$, the second integral is expanded using Eq. (B3). For the first integral in ${\mathrm{lte}}_C^{\text{Rev}}$, we use Eq. (A1) with F=C, q=q_C, and $\mathit{\varphi }=q\left(t_n\right)+\mathrm{\Delta }t(A^{*}+B^{*})$ so that $F\left(q\right[\mathit{\eta };\mathit{\varphi }\left]\right)=C\left(q_C\right[\mathit{\eta };q\left(t_n\right)+\mathrm{\Delta }t(A^{*}+B^{*})\left]\right)$, to find

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& C\left(q_C\right[\mathit{\eta };q\left(t_n\right)+\mathrm{\Delta }t(A^{*}+B^{*})\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tC\left(q\right(t_n)+\mathrm{\Delta }t(A^{*}+B^{*}\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}\left(q\right(t_n)+\mathrm{\Delta }t(A^{*}+B^{*}\left)\right)\\ & \times {\displaystyle \frac{\mathrm{d}{q}_C}{\mathrm{d}t}}[\mathrm{0};q(t_n)+\mathrm{\Delta }t(A^{*}+B^{*}\left)\right]+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),\end{array}$$

which can be simplified using Eq. (6) to get

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& C\left(q_C\right[\mathit{\eta };q\left(t_n\right)+\mathrm{\Delta }t(A^{*}+B^{*})\left]\right)\mathrm{d}\mathit{\eta }=\\ & \mathrm{\Delta }tC\left(q\right(t_n)+\mathrm{\Delta }t(A^{*}+B^{*}\left)\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}C\right){|}_{q=q\left(t_n\right)+\mathrm{\Delta }t(A^{*}+B^{*})}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}\end{array}$$

To continue the expansion, use Eqs. (B4) and (B5) to see that

$$\begin{array}{rl}q\left(t_n\right)& +\mathrm{\Delta }t(A^{*}+B^{*})=q\left(t_n\right)+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }\\ & +\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }\\ & =q\left(t_n\right)+\mathrm{\Delta }t\left(A\right(q\left(t_n\right))+B(q\left(t_n\right)\left)\right)\\ & +\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{2}}),\end{array}$$

which gives

$$\begin{array}{rl}\mathrm{\Delta }tC& \left(q\right(t_n)+\mathrm{\Delta }t(A^{*}+B^{*}\left)\right)=\mathrm{\Delta }tC\left(q\right(t_n\left)\right)\\ & +(\mathrm{\Delta }t{)}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}(A+B)\right){|}_{q=q\left(t_n\right)}+\mathcal{O}(\left(\mathrm{\Delta }t{)}^{\mathrm{3}}\right)\end{array}$$

and

$$\begin{array}{rl}{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}C\right){|}_{q=q\left(t_n\right)+\mathrm{\Delta }t(A^{*}+B^{*})}& ={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}C\right){|}_{q=q\left(t_n\right)}\\ & +\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}$$

This allows the simplification of Eq. (B6) to

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& C\left(q_C\right[\mathit{\eta };q\left(t_n\right)+\mathrm{\Delta }t(A^{*}+B^{*})\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tC\left(q\right(t_n\left)\right)\\ & +(\mathrm{\Delta }t{)}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}(A+B)\right){|}_{q=q\left(t_n\right)}+{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}C\right){|}_{q=q\left(t_n\right)}\\ & +\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}$$

The expansions of the integrals in ${\mathrm{lte}}_C^{\text{Rev}}$ are now combined to find

$${\displaystyle }{\mathrm{lte}}_C^{\text{Rev}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}(A+B)\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),$$

which is equivalent to the expression in Eq. (31).

B2 Original splitting

The mapping in the original splitting described in Eq. (20) can be written as

$$\begin{array}{rl}{\mathcal{F}}_{\mathrm{\Delta }t}^{\text{Ori}}& \left(q\right(t_n\left)\right)=q_C[\mathrm{\Delta }t;q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right]\\ & =q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]\\ & +\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}C\left(q_C\right[\mathit{\eta };q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]\left]\right)\mathrm{d}\mathit{\eta }\\ & =q_A[\mathrm{\Delta }t;q(t_n\left)\right]+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\mathrm{d}\mathit{\eta }\\ & +\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}C\left(q_C\right[\mathit{\eta };q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]\left]\right)\mathrm{d}\mathit{\eta }\\ & =q\left(t_n\right)+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }\\ & +\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\mathrm{d}\mathit{\eta }\\ & +\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}C\left(q_C\right[\mathit{\eta };q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]\left]\right)\mathrm{d}\mathit{\eta }.\end{array}$$

Thus, the local truncation error is expressed as

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}{\mathcal{F}}_{\mathrm{\Delta }t}^{\text{Ori}}& \left(q\right(t_n\left)\right)-q\left(t_{n+\mathrm{1}}\right)\mathrm{d}\mathit{\eta }=\underset{{\mathrm{lte}}_A^{\text{Ori}}}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q_A\right[\mathit{\eta };q\left(t_n\right)\left]\right)-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}A\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }}}\\ & +\underset{{\mathrm{lte}}_B^{\text{Ori}}}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }}}\\ & +\underset{{\mathrm{lte}}_C^{\text{Ori}}}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}C\left(q_C\right[\mathit{\eta };q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]\left]\right)\mathrm{d}\mathit{\eta }-\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}C\left(q\right[\mathit{\eta };q\left(t_n\right)\left]\right)\mathrm{d}\mathit{\eta }}}.\end{array}\end{array}$$

Note that ${\mathrm{lte}}_A^{\text{Ori}}$ is equivalent to ${\mathrm{lte}}_A^{\text{Rev}}$, which has already been derived in Sect. B1 and is equivalent to the expression in Eq. (25). For ${\mathrm{lte}}_B^{\text{Ori}}$, the second integral is expanded using Eq. (B2). For the first integral in ${\mathrm{lte}}_B^{\text{Ori}}$, use Eq. (A1), with F=B, q=q_B, and $\mathit{\varphi }=q_A[\mathrm{\Delta }t;q(t_n\left)\right]$ so that $F\left(q\right[\mathit{\eta };\mathit{\varphi }\left]\right)=B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]$, to find

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}B& \left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tB\left(q_A\right[\mathrm{\Delta }t;q\left(t_n\right)\left]\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}\left(q_A\right[\mathrm{\Delta }t;q\left(t_n\right)\left]\right){\displaystyle \frac{\mathrm{d}{q}_B}{\mathrm{d}t}}[\mathrm{0};q_A[\mathrm{\Delta }t;q\left(t_n\right)]\\ & +\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),\end{array}$$

which can be simplified using Eq. (6) to get

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tB\left(q_A\right[\mathrm{\Delta }t;q\left(t_n\right)\left]\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}B\right){|}_{q=q_A[\mathrm{\Delta }t;q(t_n\left)\right]}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}$$

We can now use the expansion of q_A[Δt;q(t_n)] in Eq. (A9) to simplify further:

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& B\left(q_B\right[\mathit{\eta };q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tB\left(q\right(t_n\left)\right)\\ & +(\mathrm{\Delta }t{)}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}A\right)+{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}B\right){|}_{q=q\left(t_n\right)}+\mathcal{O}(\left(\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

The expansions of the integral in ${\mathrm{lte}}_B^{\text{Ori}}$ are now combined to find

$${\displaystyle }{\mathrm{lte}}_B^{\text{Ori}}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}B}{\mathrm{d}q}}(A-C)\right){|}_{q=q\left(t_n\right)}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),$$

which is equivalent to the expression in Eq. (26). For ${\mathrm{lte}}_C^{\text{Ori}}$, the second integral is expanded using Eq. (B3). For the first integral in ${\mathrm{lte}}_C^{\text{Ori}}$, use Eq. (A1), with F=C, q=q_C, and $\mathit{\varphi }=q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]$ so that $F\left(q\right[\mathit{\eta };\mathit{\varphi }\left]\right)=C\left(q_C\right[\mathit{\eta };q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]\left]\right)$, to find

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& C\left(q_C\right[\mathit{\eta };q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tC\left(q_B\right[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}\left(q_B\right[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right){\displaystyle \frac{\mathrm{d}{q}_C}{\mathrm{d}t}}[\mathrm{0};q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right]\\ & +\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}),\end{array}$$

which can be simplified using Eq. (6) to get

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& C\left(q_C\right[\mathit{\eta };q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]\left]\right)\mathrm{d}\mathit{\eta }=\mathrm{\Delta }tC\left(q_B\right[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)\\ & +{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}C\right){|}_{q=q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}\end{array}$$

To continue the expansion, use

$$\begin{array}{rl}{q}_B& [\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]=q_B[\mathrm{0};q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]\\ & +\mathrm{\Delta }t{\displaystyle \frac{\mathrm{d}{q}_B}{\mathrm{d}t}}[\mathrm{0};q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{2}})\\ & =q_A[\mathrm{\Delta }t;q(t_n\left)\right]+\mathrm{\Delta }tB\left(q_A\right[\mathrm{\Delta }t;q\left(t_n\right)\left]\right)+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{2}})\end{array}$$

to get

$$\begin{array}{rl}\mathrm{\Delta }t& C\left(q_B\right[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)=\mathrm{\Delta }tC\left(q_A\right[\mathrm{\Delta }t;q\left(t_n\right)\left]\right)\\ & +\left(\mathrm{\Delta }{t}^{\mathrm{2}}\right)\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}B\right){|}_{q=q_A[\mathrm{\Delta }t;q(t_n\left)\right]}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}})\end{array}$$

and

$$\begin{array}{rl}{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}& \left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}C\right){|}_{q=q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]}\\ & ={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}C\right){|}_{q=q_A[\mathrm{\Delta }t;q(t_n\left)\right]}+\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}$$

We can use the expansion of q_A[Δt;q(t_n)] in Eq. (A9) to further simplify the above terms to the following:

$$\begin{array}{rl}\mathrm{\Delta }t& C\left(q_B\right[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q(t_n\left)\right]\left]\right)=\mathrm{\Delta }tC\left(q\right(t_n\left)\right)\\ & +(\mathrm{\Delta }t{)}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}A\right){|}_{q=q\left(t_n\right)}+(\mathrm{\Delta }{t}^{\mathrm{2}})\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}B\right){|}_{q=q\left(t_n\right)}+\mathcal{O}(\left(\mathrm{\Delta }t{)}^{\mathrm{3}}\right)\\ & =\mathrm{\Delta }tC\left(q\right(t_n\left)\right)+(\mathrm{\Delta }t{)}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}(A+B)\right){|}_{q=q\left(t_n\right)}+\mathcal{O}(\left(\mathrm{\Delta }t{)}^{\mathrm{3}}\right)\end{array}$$

and

$$\begin{array}{rl}{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}& \left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}C\right){|}_{q=q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]}={\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}C\right){|}_{q=q\left(t_n\right)}\\ & +\mathcal{O}\left(\right(\mathrm{\Delta }t{)}^{\mathrm{3}}).\end{array}$$

This allows the simplification of Eq. (B8) to

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}& C\left(q_C\right[\mathit{\eta };q_B[\mathrm{\Delta }t;q_A[\mathrm{\Delta }t;q\left(t_n\right)\left]\right]\left]\right)\mathrm{d}\mathit{\eta }\\ & =\mathrm{\Delta }tC\left(q\right(t_n\left)\right)+{\displaystyle \frac{(\mathrm{\Delta }t{)}^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}C\right){|}_{q=q\left(t_n\right)}\\ & +(\mathrm{\Delta }t{)}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}(A+B)\right){|}_{q=q\left(t_n\right)}+\mathcal{O}(\left(\mathrm{\Delta }t{)}^{\mathrm{3}}\right).\end{array}$$

The expansions of the integrals in ${\mathrm{lte}}_C^{\text{Ori}}$ are now combined to find

$${\displaystyle }{\mathrm{lte}}_C^{\text{Ori}}=(\mathrm{\Delta }t{)}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{d}C}{\mathrm{d}q}}(A+B)\right){|}_{q=q\left(t_n\right)}+\mathcal{O}(\left(\mathrm{\Delta }t{)}^{\mathrm{3}}\right),$$

which is equivalent to the expression in Eq. (27).