Interactive comment on “ Numerical issues associated with compensating and competing processes in climate models : an example from ECHAM-HAM ”

1) The description of scheme HAM1 should be more precise, it is unclear how the limiter is exactly implemented and whether the clipping acts on S** or on S_t+deltat or on both. Furthermore, it is not clearly specified if the 95% limiting is also used in scheme 1EP. In general, it would be better to express all these limiting steps by formulae, including appropriate operators like max(S,0) in the definition of the numerical method.


Introduction
In the past decades, the climate modeling community has been moving gradually towards high-resolution and processbased modeling.More and more complex processes such as the details of aerosol lifecycle and cloud microphysics are be-ing brought into global and regional climate models.During this evolution, the range of timescales explicitly represented in the models is broadening substantially.From a mathematical point of view, this means the system of differential equations is not only expanding but in the meanwhile getting considerably stiffer, posing great challenges to the numerical methods applied in climate models.
Traditionally, numerics is considered by many as the main focus of dynamical core developers and not so much for physicists who design parameterization schemes for sub-grid processes.The air pollution and chemistry transport modeling community, as well as various numerical weather forecast centers, have paid substantial attention to the use of numerical techniques in complex models (e.g., Girard and Delage, 1990;Beljaars, 1991;Teixeira, 2000;Benard et al., 2000;Jacobson, 2002;Beljaars et al., 2004;Wood et al., 2007;Zaveri et al., 2008;Schlegel et al., 2012;Tudor, 2012), while climate modelers have focused less on this issue.In this paper we present an example from the global aerosol-climate model ECHAM-HAM (Stier et al., 2005;Zhang et al., 2012) to demonstrate that numerical errors associated with stiff systems can lead to significant systematic biases in simulations at typical spatial and temporal scales considered in climate research.The example is also relevant to the numerical treatment of many other processes in climate models.
What motivated this work was the dramatic changes in the sulfuric acid (H 2 SO 4 ) gas concentration when ECHAM-HAM was upgraded from version 1 to version 2 (hereafter referred to as HAM1 and HAM2, respectively).The use of a new time stepping scheme for the H 2 SO 4 gas equation resulted in a significant increase of H 2 SO 4 gas concentration Published by Copernicus Publications on behalf of the European Geosciences Union.
at most grid points (Zhang et al., 2012).On the one hand, the new scheme outperforms the old one in box model tests presented by Kokkola et al. (2009).On the other hand, comparisons of the global model simulation with (the still rare) H 2 SO 4 gas measurements seem to suggest that the new model version is associated with larger positive biases (D.O'Donnell, personal communication, 2012).In this paper we carry out convergence tests for the two time stepping schemes and analyze the characteristics of the numerical errors.The aim is to identify a better scheme from a numerical perspective, and quantify the remaining biases.Impact of these biases on the simulated aerosol and cloud condensation nuclei (CCN) number concentration is also discussed.
As is elaborated later in the paper, the key to an accurate numerical solution of the H 2 SO 4 gas equation is the proper balances between strongly compensating and competing processes.Such situations of process interaction are often encountered in other components of the climate models as well.Examples include the liquid water budget in cloud microphysics (P. Caldwell, personal communication, 2013) and the wind profile in the near-surface levels (Beljaars, 1991).In this paper we consider the sulfuric acid gas equation as a prototype problem of this kind.By testing several time stepping methods in addition to those used in HAM1 and HAM2, we attempt to obtain some generally useful conclusions regarding the numerical representation of process interactions in climate models.
The remainder of the paper is organized as follows: Sect. 2 introduces the sulfuric acid gas equation in ECHAM-HAM and outlines the time stepping methods.Section 3 presents results of the convergence test and establishes the reference solution.Section 4 focuses on the issue of strongly compensating processes, and Sect. 5 the competing processes.Conclusions from this study are summarized in Sect.6.

Methodology
This section first briefly introduces the ECHAM-HAM model to set the context.The sulfuric acid gas equation is then described, together with the time stepping schemes used in HAM1 and HAM2.Other integration schemes used in the sensitivity experiments are also introduced.Thereafter, the simulations for testing these schemes are described.

Overview of ECHAM-HAM
ECHAM-HAM is a global aerosol-climate model developed for understanding the distribution, properties and lifecycle of tropospheric aerosols as well as their interactions with climate.
The atmospheric general circulation model ECHAM5 (Roeckner et al., 2003(Roeckner et al., , 2006) ) solves the hydrostatic equations of fluid motion using the spectral transform method with triangular truncation.In the vertical, the model uses a pressure-based terrain-following coordinate with discretization methods following Simmons and Burridge (1981).The highest computational level is located at 10 hPa.The large-scale transport of water substances and other tracers is handled by the flux-form finite volume algorithm of Lin and Rood (1996), assuming the fields vary with parabolic sub-grid distributions.Cumulus convection and convective tracer transport are described by the mass flux scheme of Tiedtke (1989), with further modifications by Nordeng (1994).The turbulent transport of momentum, heat, moisture and tracers is represented by the eddy-diffusivity scheme of Brinkop and Roeckner (1995) which involves a prognostic equation for the turbulent kinetic energy.Short-wave and long-wave radiative transfer calculations follow the methods of Fouquart and Bonnel (1980) and Mlawer et al. (1997), respectively.
The aerosol module HAM was first developed by Stier et al. (2005), and has gone through various updates in recent years.A summary can be found in Zhang et al. (2012).The aerosol population in the atmosphere is described by the mass concentrations of different chemical species (sulfate, black carbon, organic matter, sea salt, and mineral dust), and the number concentrations of different particle size classes.The particle size distribution is mathematically described by 7 log-normal modes, 4 of which correspond to soluble aerosols that are internally mixed (meaning one particle can contain more than one chemical composition), while the other 3 modes are insoluble, and externally mixed (i.e., each particle contains only one chemical species).In the model, aerosols can form in the atmosphere through nucleation processes.The mechanisms considered in HAM2 include neutral and charged nucleation of H 2 SO 4 and H 2 O (Kazil and Lovejoy, 2007;Kazil et al., 2010).In the planetary boundary layer over forested areas, the nucleation of H 2 SO 4 and an organic compound can be simulated by the kinetic nucleation parameterization of Laakso et al. (2004) and Kuang et al. (2008), or the cluster activation scheme of Kulmala et al. (2006) andRiipinen et al. (2007).Aerosol particles can also be directly released into the atmosphere through natural and/or anthropogenic emission.The emission fluxes are interactively calculated for sea salt and dust (Monahan et al., 1986;Smith and Harrison, 1998;Tegen et al., 2002;Cheng et al., 2008), and prescribed for the other species.Microphysical processes, such as coagulation and condensation of H 2 SO 4 gas, are considered following Vignati et al. (2004).The parameterization of aerosol water uptake is based on the work of Petters and Kreidenweis (2007).Aerosols are removed from the atmosphere by gravitational settlement (Stier et al., 2005), turbulent dry deposition (Kerkweg et al., 2006), as well as in-cloud and below-cloud scavenging (Stier et al., 2005;Zhang et al., 2012, and references therein).
Regarding the connection between aerosols and climate in ECHAM-HAM, the atmosphere model provides meteorological conditions that are needed for the aerosol calculations (emissions, microphysics, removal processes), and handles Fig. 1.Annually and zonally averaged sulfuric acid gas concentrations (unit: number of molecules per cm 3 ) simulated by two configurations of the aerosol-climate model ECHAM-HAM, version 2. Panel (a) corresponds to time stepping scheme 1 in Table 1 which solves the sulfuric acid equation (Eq. 1) using the Euler forward scheme with sequential splitting; Panel (b) uses scheme 2 in Table 1, originally introduced by Kokkola et al. (2009).The simulation setups are described in Sect.2.3.the large-scale and sub-grid-scale transport of aerosols and their precursors.Aerosols affect atmospheric circulation by changing the radiation budget and through their impacts on cloud microphysics (Lohmann et al., 2007;Lohmann and Hoose, 2009).

Sulfuric acid gas equation
In ECHAM-HAM, an ordinary differential equation (ODE) of the form is included to represent the link between sulfur chemistry and aerosol microphysics.Here S denotes the concentration of H 2 SO 4 gas in the unit of molecules per cubic centimeter.P is the source term related to chemical production and transport processes.C•S describes the condensation of H 2 SO 4 gas on pre-existing aerosol particles.N(S) denotes the H 2 SO 4 gas loss rate due to aerosol nucleation, generally a nonlinear function of S. Within each step of the time integration, the source term P and the condensation coefficient C are kept constant.
In the old model HAM1, Eq. ( 1) is solved by the Euler forward scheme using sequential splitting.The concentration S is updated in three consecutive steps, each considering one term on the r.h.s. of Eq. (1): (2) Note that in Eq. ( 3) the H 2 SO 4 gas condensation is limited to 95 % of the available amount, while in Eq. ( 4) there is a limiter applied to nucleation to avoid negative concentrations.This algorithm is referred to as "scheme 1" in the remainder of the paper.
In HAM2, a two-step time integration scheme proposed by Kokkola et al. (2009) is employed.First, production and condensation are considered together using This expression is an analytical solution of the productioncondensation equation (i.e., Eq. 1 with N (S) = 0) when assuming the production rate P and condensation coefficient C are constant within the time interval t.Subsequently, the nucleation sink N is computed using the intermediate concentration S * , and adjusted with an Euler-backward factor 1/(1 + C t). (This adjustment factor comes from an attempt to discretize Eq. ( 1) with the Euler-backward scheme in which iterative evaluations of aerosol nucleation are employed to avoid a nonlinear solver.A detailed explanation can be found in Sect. 3 of Kokkola et al. (2009).)The new concentration S t+ t is given by This integration method is referred to as "scheme 2" in the following.As can be seen in Fig. 1, the use of Eqs. ( 5) and (6) instead of the old scheme results in considerably higher H 2 SO 4 gas concentrations at most grid points in the model domain.(Details of the simulation setup are given in the next subsection.) In order to explain this difference and to analyze the properties of the two schemes, the following time stepping schemes are also tested: -Scheme 1EP: Similar to scheme 1 but using parallel splitting for production and condensation, i.e., As in scheme 1, sequential splitting is used between production/condensation and nucleation.Equations ( 7) www.geosci-model-dev.net/6/861/2013/Geosci.Model Dev., 6, 861-874, 2013 and ( 8) also include the same limiters as in scheme 1 (Eqs.2-4), for the purpose of a clean comparison between the sequential and parallel splitting of production and condensation.
-Scheme 1Im: Similar to scheme 1EP but replace Eq. ( 7) by the trapezoidal implicit scheme A limiter S * = max(0, S * ) is applied, assuming that all available H 2 SO 4 gas, rather than 95 % of it, can condense on existing aerosol particles.As discussed later in the last paragraph of Sect.4.2, this change in the limiter helps to eliminate positive errors in regions of high aerosol loading.
-Scheme 2CP: Use analytical solution for the production-condensation equation as in schemes 2 and 2C (Eq.5), but parallel splitting between production/condensation and nucleation, i.e., Eq. ( 5) followed by The time integration methods described above, except for scheme 2CP, are based on sequential splitting between nucleation and the rest of the ODE, limiting the numerical convergence to first order.In addition to these schemes, we evaluate two methods that solve production, condensation and nucleation simultaneously.These schemes are based on a Taylor expansion of the nucleation sink Substituting Eq. ( 11) into (1), we get a linearized differential equation with Now that Eq. ( 12) has constant coefficients within one time step, we can apply the analytical solution to get The derivation of Eq. ( 12) can be interpreted as linearizing the right-hand side of Eq. ( 1) around the initial condition of the time step.This is one of the the essential ideas behind the widely used Rosenbrock methods (Rosenbrock, 1963;Hairer and Wanner, 2004).It can be shown that Eq. ( 15) is equivalent to the so-called exponential Rosenbrock-Euler method (Hochbruck et al., 2009;Schweitzer, 2013) which provides second-order accuracy.We refer to Eq. ( 15) as "scheme 3A".
For comparison with the first-order schemes outlined earlier, we also test a "scheme 3B" that solves Eq. ( 12) using the Euler-backward method, i.e., with the limiter It can be shown that Eq. ( 16) is equivalent to the one-stage Rosenbrock method (Rosenbrock, 1963;Hairer and Wanner, 2004).
Because the derivative (dN/dS) t is usually not readily provided by the parameterization, and often nontrivial to derive from the original formulation unless the scheme is relatively simple, it needs to be estimated numerically.We have tested the approximation with various values for β, and found the results to be rather insensitive.Therefore, we present here the simulations performed with β = 0 which requires only one calculation of the nucleation sink per time step.The use of β = 0 for Eq. ( 18) and the Euler-backward scheme (16) for Eq. ( 12) effectively gives the solver of Jacobson (2002) (Eq.( 16) and paragraph 25 therein; see also Eqs. (16.68) and (16.74) in Jacobson, 2005).As pointed out by Jacobson (2002) and shown later in Sect. 5 of this paper, solving nucleation and condensation together helps to correctly represent the competition between the two processes for the available sulfuric acid gas.In all numerical schemes described above, the production rate P and condensation coefficient C in the H 2 SO 4 gas Eq. ( 1) are frozen within one model time step.It is believed that for global model simulations, this assumption is justified when the physics time steps is on the order of one hour or shorter, for the following reasons: the condensation coefficient C is mainly determined by the surface area of pre-existing aerosols, which in turn is a function of aerosol concentration and size distribution.The typical aerosol lifetime is about 0.6 days for sea salt and 4-7 days for other species, much longer than the lifetime of H 2 SO 4 gas (0.01-0.02 days).The chemical production rate P is primarily determined by the precursor and oxidant (i.e., SO 2 and OH) concentrations and the ambient temperature.SO 2 gas has a typical lifetime of 1-2 days, while for OH the ECHAM-HAM model considers only seasonal cycle and diurnal cycle.Therefore, both P and C are expected to vary relatively smoothly with time in comparison to the H 2 SO 4 gas itself, justifying the use of a frozen coefficient in Eq. ( 1).In principle one could drop the assumption of constant C and solve Eq. ( 1) in combination with the aerosol equations.But that would lead to a very complicated system, given the large number of prognostic variables (for the mass and number concentrations of different aerosol species and size ranges) and the numerous parameterized, highly nonlinear microand macro-physical processes involved.To the best of our knowledge, it is common for aerosol-climate models and aerosol-chemistry models to solve gas condensation equations assuming constant coefficients within one time step (e.g., Jacobson, 2002;Zaveri et al., 2008), as we do here.Considering that the surface areas of small particles are more readily affected by gas condensation and aerosol nucleation in comparison with large particles, one could alternatively couple Eq. ( 1) with concentration equations of the nucleation mode aerosols, but consider only the changes in aerosol mass due to gas condensation and change in aerosol number due to nucleation.Such a equation system may give more accurate results in areas where the nucleation mode particles dominate the aerosol population (e.g., in the tropical tropopause).This alternative is worth evaluating in the future but not investigated here.In this study we focus on solving Eq. ( 1) in an isolated manner, and try to answer the question "given Eq. ( 1) with constant P and C, how does process splitting affect the solution of the equation and the aerosol concentration in the model".
The test strategy employed in this paper is to carry out simulations with ECHAM-HAM, rather than a box model, in order to evaluate the numerics in different regions and regimes in real-world simulations.Comparison with observation is not presented in this paper because we want to focus on numerical issues and separate them from the influence of parameterization schemes as well as model biases from other sources.The reference solution is established by carrying out convergence tests using small time steps for the sulfuric acid equation.Time steps of the rest of the model stay unchanged.

Simulations
Global simulations with ECHAM-HAM2 are performed for the year 2000.The model system is forced by the sea surface temperature distribution and sea ice concentrations compiled by the Second Atmospheric Model Intercomparison Project (http://www-pcmdi.llnl.gov/projects/amip/). Aerosol emissions are specified according to Dentener et al. (2006), except that the formation of secondary organic aerosol is explicitly represented (O'Donnell et al., 2011).The model meteorology is nudged towards the ERA-40 reanalysis (Uppala et al., 2005).These set-ups are essentially the same as in Zhang et al. (2012).In our simulations, the boundary layer aerosol nucleation scheme of Laakso et al. (2004) and Kuang et al. (2008) is switched on.
Although HAM2 is most often run at T63L31 resolution (approximately 2 • grid spacing in the horizontal, 31 vertical levels), we use T42L19 (approximately 3 • , with 19 vertical levels) in this work, due to the large number of simulations involved.A subset of the experiments has been repeated at T63L31, in which it was found that although the absolute values of the numerical error are generally smaller than those at T42L19 (as expected), the relative differences between results from different numerical schemes are similar to what is presented here.The main conclusions from our investigation are not affected by the choice of model resolution.
For clarity we note that the dynamical core of ECHAM uses a leap-frog integration method with semi-implicit treatment for linear gravity waves.The default time step is 30 min at resolution T42L19.The physics parameterizations use two-time-level schemes that advance the model state from t − t D to t + t D where t D stands for the time step of the dynamical core.The t used in Eqs. ( 2)-( 16) is thus equal to 2 t D .

Convergence test and reference solution
Convergence tests are performed for the time integration schemes described in Sect.2.2 using up to 256 sub-steps per each physics step.The resulting annual mean H 2 SO 4 gas burden is presented in Table 1.It can be seen that discrepancies caused by the use of different time stepping schemes decrease when more sub-steps are used.With 128 and 256 sub-steps (28 s and 14 s sub-step size, respectively), results from the seven simulations 1EP, 1Im, 2, 2C, 2CP, 3A and 3B are less than 0.2 % apart.Results given by the HAM1 scheme are less than 2 % (with 128 sub-steps) and 1 % (with 256 sub-steps) different from the other simulations.Based on this table, we choose to use the average of the simulations in the rightmost column (excluding scheme 1) as the reference solution in further analysis.
Relative differences with regard to the reference solution in the globally integrated annual mean H 2 SO 4 gas concentration, condensation rate and nucleation sink in the simulations listed in Table 1 are shown in Fig. 2. In addition to a clear trend of convergence as the number of sub-steps increases, it can be seen that scheme 2 (light green line) is associated with considerably smaller errors than scheme 1 (dark blue line) in all three quantities when the number of sub-steps is small.We thus conclude that the change of time stepping scheme from HAM1 to HAM2 is a solid improvement, as least from the point of view of numerics.Regarding the concern with possibly larger positive biases in H 2 SO 4 gas concentration in HAM2, Zhang et al. (2010) noticed that the SO 2 oxidation in HAM was considerably stronger than in a different aerosol model.Figure 2b   The HAM2 time integration scheme (Kokkola et al., 2009).( 5), ( 6) 1.3848 1.3582 1.3344 1.3346 2C The HAM2 scheme but without the Euler backward correction of nucleation.
(5), ( 8 gas loading down to a satisfactory level.In the new model version, such cancellation between physical and numerical errors no longer exists.In order to constrain the H 2 SO 4 gas loading in HAM2, it will probably be useful to re-evaluate the H 2 SO 4 gas production and related parameterizations in ECHAM-HAM. In the following sections, further analyses of the simulations shown in Table 1 and Fig. 2 are presented.Before going into the details, we show in Fig. 3 the zonal and annual mean vertical cross sections of the H 2 SO 4 gas concentration as well as its source and sinks.The main locations of sulfate particle formation are the tropical tropopause (Fig. 3b), where high relative humidity and low temperature provide favorable conditions for the neutral nucleation of H 2 SO 4 and H 2 O, and the near-surface layers (Fig. 3b) where the concentrations of H 2 SO 4 gas are high.Chemical production and condensation of the H 2 SO 4 gas peak in the near-surface levels where both the precursors and aerosol particles are abundant.
A comparison of the magnitude of the source and sink terms in Fig. 3b-d reveals that production and condensation almost compensate each other.In terms of the long-term H 2 SO 4 gas budget, the nucleation sink is effectively balanced by the residual of two large terms of opposite signs.This is a typical case in numerical modeling where computational error easily contaminates the results.Further discussions on this issue are included in Sect. 4. In the tropical tropopause, nucleation plays a non-negligible role in the H 2 SO 4 gas budget where its magnitude can exceed 20 % of that of condensation in terms of annual and zonal mean (Fig. 3e).In this region, the partitioning of available sulfuric acid gas between condensation and nucleation, or in other words, the competition between the two processes, becomes important for an accurate representation of nucleation.This is further discussed in Sect. 5.

Operator splitting
Based on the previous section, we can now explain the increase of H 2 SO 4 gas concentration from Fig. 1a, b.The key reason is the replacement of a sequential splitting between production and condensation by the analytical solution (Eq. 5) of the production-condensation equation.
In the real atmosphere, production and condensation occur simultaneously.Condensation prevents H 2 SO 4 gas from linearly increasing, and forces it to asymptotically approach the equilibrium concentration P /C with an e-folding time of C −1 (cf.Eq. 5).A sequential splitting scheme first updates H 2 SO 4 gas by considering only the production, resulting in a positive error in the intermediate concentration that is used for computing the condensation rate.The scheme thus features systematic overestimate of condensation and negative bias in H 2 SO 4 gas burden, as can been seen in Fig. 2a, b (dark blue lines).The discretization errors are particularly large when the integration time step is long in comparison to the condensation time scale.
To verify this reasoning, we carried out two simulations, 1EP and 1Im, and compare them with simulation 1.The explicit parallel splitting of production and condensation is based on the thinking that since the two terms are largely compensating each other, the gas concentration will not change dramatically in one time step, thus the step-average concentration can be approximated by the initial value (Euler forward scheme, first-order accuracy).Parallel splitting reduces the absolute error in H 2 SO 4 gas burden by a factor of 10 (Fig. 2a, purple vs. dark blue line).However, the condensation rate features a considerable negative bias (Fig. 2b, purple line) because the condensation of newly produced H 2 SO 4 gas is not considered.The implicit scheme, in contrast, considers production and condensation together, and approximates the step-average concentration by the arithmetic average of the initial and ending values.This method turns out much more accurate than schemes 1 and 1EP.The globally and annually averaged H 2 SO 4 gas concentration and condensation rate are very close to those obtained with the analytical solution (Fig. 2a, b, scheme 1Im vs. 2C).These results suggest that for strongly compensating processes, both sequential parallel splitting can lead to large numerical errors when used in combination with long time steps.Although our simulations 1 and 1EP both use Euler forward time stepping, the explicit nature of the schemes is not the main cause of errors in this case.In the aerosol module MAM3 of the CAM5 model (Liu et al., 2012), the production-condensation equation is solved with a sequential splitting method like in HAM1, but using the analytical solution of the condensation equation dS/dt = −C • S, i.e., Eq. ( 3) is replaced by Despite the fact that Eqs. ( 2) and ( 19) are exact solutions of the production and condensation equations, respectively, results from the MAM3 module also have severe negative biases in H 2 SO 4 gas concentration in the near-surface levels (Liu and Easter, personal communication, 2012).The real culprit of the systematic errors in our simulations 1, 1EP and in MAM3 is the use of operator splitting with a large time step.The HAM2 method and the implicit scheme 1Im solve production and condensation together, thus produce more accurate results.

Comments on sub-stepping and clipping
Many of the parameterization schemes in contemporary climate models are nonlinear, which can make it impractical or inconvenient to use analytical solutions and/or implicit methods.In such cases, sub-stepping is often used to handle fast processes.Here we want to point out several caveats related to the use of sub-stepping.In Fig. 2a, the H 2 SO 4 gas burden errors obtained with scheme 1EP are of opposite signs when 1 and 4 sub-steps are used.This reflects the inherent nature of the explicit scheme.According to the stability analysis in, e.g., Chapter 2 of Butcher (2008), the Euler forward scheme leads to oscillatory behavior in the solution of the production-condensation equation when the step size exceeds the e-folding time 1/C,  (b-d) Relative differences with regard to the reference solution in simulations using scheme 1 (cf.Sect.2.2), and using scheme 1EP without and with 2 sub-steps.All panels are plotted for the lowest model level.Similar figures of the other months convey the same message, although the spatial distributions are different because of the seasonal cycle.and causes instability when t > 2/C.In Fig. 4 the Januarymean near-surface condensation coefficients are shown for the reference solution established in Sect.3, which indicates that the 60 min step size (default for the parameterized physics at T42L19) is unstable for most of the grid points in the lower troposphere.Such instability does not cause floating point exception in our simulations, or H 2 SO 4 gas concentrations that are orders of magnitude off, thus can be overlooked.The resulting numerical errors, however, are substantial.Similar instability problems also exist for other explicit Fig. 6.Relative errors of H 2 SO 4 gas concentration in box model calculations that solve the production-condensation equation using the Euler forward scheme for an 1 h integration time.Typical values of the production rate P , condensation coefficient C, and initial value S 0 in Southeast China are used, which represents region associated with persistent positive errors in Fig. 5b-e.The dashed green line indicate the stability threshold t = 2/C (cf., e.g., Chapter 2 in Butcher, 2008).In panel (a) the estimated condensation rate at each step of time integration is limited to 95 % of the total available H 2 SO 4 gas, while in panel (b) this clipping factor is set to 1.Note that the two panels use different scales for the y-axis.Further details can be found in Sect.4.2.methods, e.g., Runge-Kutta and explicit predictor-corrector methods.Sensitivity experiments have been performed but are not shown here.
To reduce the errors, sub-stepping is needed to ensure sufficiently small step size.In current climate models, this is commonly employed with a fixed number of steps determined by investigations of simplified test cases or evaluation of zonal mean statistics (e.g., Posselt and Lohmann, 2008;Morrison and Gettelman, 2008;Gettelman et al., 2008).In our simulations with the 1EP scheme, 8 sub-steps (7.5 min sub-step size) turns out sufficient to provide a less than 2 % error in the annual mean H 2 SO 4 gas burden, and less than 15 % errors in the annually averaged zonal mean concentration.However, the instantaneous values of condensation coefficients can be so high as to require more than 200 sub-steps.This is problematic for studies that use a fixed small number of sub-steps in applications that investigate re-gional features and impacts.Given that very small step sizes (e.g., on the order of seconds) are expensive to use globally and also unnecessary for grid points with relatively weak condensation, the use of dynamically controlled time steps is worth considering.Zaveri et al. (2008) developed an adaptive time stepping scheme using a priori estimates of step size, while Herzog et al. ( 2004) used time steps dynamically adjusted according to an a posteriori error estimate.
For the H 2 SO 4 gas equation in ECHAM-HAM, we tested a version of the 1EP scheme in which the sub-step size is calculated from the condensation coefficient.The number of sub-steps is chosen such that the stability factor C t does not exceed unity.In order to stay with the original coding structure in ECHAM-HAM, the same number of sub-steps is applied to all grid boxes in the same CPU, determined by the largest condensation coefficient among these grid boxes.Results show that for a one-year simulation, the globally averaged number of sub-steps is about 20 when using 32 CPUs.The simulated annual mean H 2 SO 4 gas burden and condensation rate are less than 1 % different from the reference solution.
As a side remark, we note that the use of adaptive substepping can cause load balancing issues in parallel computing, when the same number of grid boxes are assigned to each CPU.This was not a problem in our simulations because the H 2 SO 4 gas equation only took a very small fraction of the total computing time.The boundary layer nucleation scheme has a rather simple formulation, and the Kazil and Lovejoy (2007) parameterization is implemented as a look-up table (Kazil et al., 2010).For more complex and computationally expensive parameterization-like cloud macro-and microphysics, the load imbalance can be significant or even become a bottle neck in parallelization.In such a case, the use of advanced domain decomposition algorithms, such as the METIS graph partitioning tool (Karypis andKumar, 1995, 1998), will probably be helpful.
Another feature worth noting in our simulations of the H 2 SO 4 gas is that despite the oscillatory convergence of the 1EP scheme and the generally negative biases in the H 2 SO 4 gas concentration simulated by scheme 1, the regions with highest aerosol loading are found to be associated with persistent positive biases in the near-surface H 2 SO 4 gas concentration, as shown in Fig. 5b-d.These large positive errors occur because as in HAM1, the condensation sink is limited to 95 % of the available H 2 SO 4 gas.This means in the case of strongly compensating production and condensation, the production is allowed to slightly exceed condensation.Sensitivity experiments indicate that if a factor of 95 % is used in scheme 1Im, similar positive errors occur (not shown).Box model calculations displayed in Fig. 6 confirm that when the clipping factor is changed from 95 % to 100 % in scheme 1EP, the persistent positive errors associated with the small sub-step numbers disappear.Meanwhile the positive errors become significantly larger (Fig. 6b vs. a) because complete depletion of sulfuric acid gas in one (sub-)step leads to a severely underestimated condensation rate in the next (sub-)step.The comparison of Fig. 6a with b suggests that the 95 % clipping factor does help to reduce errors in the heavily polluted regions when a small number of sub-steps are used.On the other hand, additional box model calculations (not shown) and the error patterns displayed in Fig. 5 indicate that the sign and magnitude of the errors depend on the characteristic condensation coefficient, thus have a strong regional variation.This again would undermine the accuracy of studies that used such a clipping when studying regional features and impacts.To reduce the impact of clipping, smaller time steps and/or more accurate time integration schemes should be applied.

Aerosol nucleation and cloud condensation nuclei
Numerical error in the simulated H 2 SO 4 gas nucleation sink comes from two sources: (i) the H 2 SO 4 gas concentration provided as input to the nucleation parameterization scheme, e.g., S * * in Eq. ( 4); and (ii) the time integration method used for the nucleation process, e.g., the correction factor 1/(1 + C t) in Eq. ( 6).The first source explains the 75 % negative bias in global mean nucleation sink given by the HAM1 numerics (scheme 1 in Fig. 2c, dark blue line), as well as the positive biases in the near-surface levels that can be seen from the zonally and annually averaged vertical cross section in Fig. 7b.The impact of source (ii) can be investigated by comparing schemes 2, 2C and 2CP in Table 1, in which the same numerical treatment is applied to production and condensation, while the coupling with nucleation is varied.The simulated zonal and annual mean nucleation sink is shown in Fig. 7c-h.
Scheme 2C applies a sequential splitting between production-condensation and nucleation (like in HAM2), but no special correction for nucleation.In regions where nucleation is non-negligible in magnitude (cf.Fig. 3e), the intermediate H 2 SO 4 gas concentration S * calculated with Eq. ( 5) (which ignores the nucleation sink) tends to have considerable positive biases, resulting in overpredicted nucleation.This is why the scheme produces positive bias in the nucleation sink near the tropopause (Fig. 7e, f).
The nucleation adjustment factor 1/ (1 + C t) that Kokkola et al. (2009) introduced to HAM2 (scheme 2, Eq. 6) helps to reduce the absolute bias in the upper troposphere, but tends to cause over-correction (negative biases in Fig. 7c, d).
From the perspective of competition between condensation and nucleation, the HAM2 scheme favors condensation because while the nucleation sink is scaled down, the condensation process does not have any competitor in the first step of the calculation.The nucleation adjustment factor causes considerable negative biases in aerosol nucleation in the lower troposphere (Fig. 7d) due to the large condensation coefficients and long time step.Close to the surface, the relative differences with regard to reference solution exceed −50 % in middle and low latitudes, and even −75 % in polluted regions.Consequently, not only the concentration of the nucleation mode aerosols are severely underestimated, but also the Aitken mode number concentration is associated with biases of about 20 %-70 % in North America, Europe and North Africa, and in East Asia.The latter can be seen in Fig. 8a, where the relative errors in the Aitken mode number concentration of the lowest model level are shown for scheme 2. The negative biases in the concentration of cloud condensation nuclei (CCN, diagnosed at 0.1 % supersaturation) generally exceed 5 % in these regions (Fig. 8b).Based on these results, we suggest removing the nucleation adjustment factor in future versions of the ECHAM-HAM model.
Without the adjustment factor (i.e., scheme 2C, Fig. 7f), the model provides good results in the near-surface levels but about 20 % relative errors in the nucleation sink near the tropical tropopause.If parallel splitting is applied between nucleation and production-condensation (scheme 2CP), errors of similar magnitude are produced although the signs are different (Fig. 7h).The time integration schemes 3A and 3B, which handle the three processes simultaneously, significantly improves the nucleation sink throughout the model domain.Despite the fact that scheme 3B has the same order of accuracy (first-order) as those using sequential splitting, it produces smaller errors (Figs.7j and 8c, d).Results from the second-order scheme 3A are similar to 3B, hence not shown.Because the purpose of solving the H 2 SO 4 gas equation in ECHAM-HAM is to represent its impact on aerosol formation and growth, based on Figs. 7 and 8 we consider schemes 3A and 3B as better choices for this model.
As pointed out earlier in Sect.3, the H 2 SO 4 gas budget in the middle and upper troposphere can be understood as the balance between nucleation (sink) and the residual of production minus condensation (source).The small errors associated with scheme 3B indicate again that strongly compensating processes need to be treated together by the numerical solver.

Conclusions
In this study we analyze various time integration methods for the sulfuric acid gas equation in the aerosol-climate model ECHAM-HAM.Tests of numerical convergence are performed using small time steps, from which a reference solution is established.It is found that the H 2 SO 4 gas concentration and sulfate particle formation rate provided by the time stepping scheme in HAM2 are associated with considerably smaller numerical errors than those from HAM1.The main reason for this improvement is the replacement of a sequential splitting between production and condensation by the analytical solution of the production-condensation equation which includes the dominant terms of the H 2 SO 4 gas budget.Although comparison with observations in other studies have provided hints that the simulated H 2 SO 4 gas concentration Fig. 7. Left column: zonally and annually averaged H 2 SO 4 gas loss rate due to aerosol nucleation (unit: molecules cm −3 s −1 ), simulated using various time stepping scheme with the default time step for model physics (i.e., no sub-stepping).Right column: the corresponding relative differences (unit: %) with regard to the reference solution shown in Fig. 3b.The relative differences are plotted only in regions where the nucleation sink is stronger than 1 molecule cm −3 s −1 .The numerical schemes are described in Sect.2.2 and Table 1. may have degraded from HAM1 to HAM2, our results indicate that the time stepping method should not be reverted.Instead, other components of the model, possibly the production of the H 2 SO 4 gas, need to be re-evaluated.
The time stepping scheme currently used in HAM2 includes a numerical adjustment for the H 2 SO 4 gas loss rate due to nucleation, which causes significant negative biases in the nucleation sink, aerosol number concentration, and CCN concentration in regions of strong condensation.We suggest this adjustment be removed in the future.Furthermore, two time stepping schemes based on the linearization of nucleation are tested in the paper.The schemes solve H 2 SO 4 gas production, condensation and nucleation simultaneously, and provide more accurate results than the current model version.
In modern climate models that include process-based physics parameterizations, the wide spectrum of timescales poses new challenges to the numerical methods.Our study demonstrates that in situations involving strongly compensating and/or competing processes, operator splitting in combination with long time step can cause severe numerical errors.The sub-stepping technique used with a fixed small number of sub-steps may give good results in terms of global or zonal mean, but can still cause large regional biases due to the extremely short characteristic timescales in certain regions.
The key conclusion from this work is that, in order to reduce the biases associated with process interactions, it is important to apply numerical methods that solve multiple processes simultaneously.Such methods, employed together with implicit time stepping, can provide high-accuracy re-sults without the need for very short time steps.As for substepping, the use of dynamically adjusted step size can be beneficial, and provide a good trade-off between numerical accuracy and computational efficiency.The potential issue of load balancing in parallel computing may be addressed by advanced domain decomposition algorithms.
suggests that the numerical error in condensation probably helped to bring the sulfuric acid www.geosci-model-dev.net/6

Fig. 3 .
Fig. 3. Zonally and annually averaged vertical cross sections of the H 2 SO 4 gas (a) concentration, (b) loss rate due to aerosol nucleation, (c) condensation rate, (d) production rate, in the reference solution established in Sect.3. Panel (e) shows the ratio between the nucleation sink and the condensation rate.

Fig. 4 .
Fig. 4.Reference solution of the January-mean H 2 SO 4 gas condensation coefficient (unit: s −1 ) in the lowest model layer.The numbers given in parentheses next to the color bar are the step sizes corresponding to the stability threshold t = 2/C of the Euler forward time integration scheme (cf., e.g., Chapter 2 inButcher, 2008).

Fig. 5 .
Fig. 5. (a)Reference solution of the January-mean near-surface H 2 SO 4 gas concentration (unit: molecules per cm 3 ).(b-d) Relative differences with regard to the reference solution in simulations using scheme 1 (cf.Sect.2.2), and using scheme 1EP without and with 2 sub-steps.All panels are plotted for the lowest model level.Similar figures of the other months convey the same message, although the spatial distributions are different because of the seasonal cycle.

Fig. 8 .
Fig.8.Relative differences with regard to to the reference solution (unit: %) of the annual mean number concentrations of the soluble Aitken mode aerosol (left column) and the cloud condensation nuclei at 0.1 % supersaturation (CCN 0.1 %, right column) in the lowest model layer.The values are plotted only in regions with Aitken mode concentration > 500 cm −3 STP and CCN > 50 cm −3 , respectively.STP stands for standard temperature and pressure (1013.25 hPa, 273.15 K).The upper row shows results obtained with time integration scheme 2. The lower row corresponds to scheme 3B.

Table 1 .
Annual mean H 2 SO 4 gas burden given as globally integrated mass of sulfur (unit: 10 −3 Tg) simulated by ECHAM-HAM at T42L19 resolution using different time stepping schemes and numbers of sub-steps for the sulfuric acid gas equation (Eq.1).