libcloudph + + 1 . 1 : aqueous phase chemistry extension of the Lagrangian cloud microphysics scheme

This paper introduces a new scheme available in the library of algorithms for representing cloud microphysics in numerical models named libcloudph++. The scheme extends the Lagrangian microphysics scheme available in libcloudph++ to the aqueous phase chemical processes occurring within cloud droplets. The representation of chemical processes focuses on the aqueous phase oxidation of the dissolved SO2 by O3 and H2O2. The Lagrangian Microphysics and Chemistry (LMC) scheme allows tracking the changes in the cloud condensation nuclei (CCN) distribution caused by both collisions between 5 cloud droplets and aqueous phase oxidation. The scheme is implemented in C++ and equipped with bindings to Python which allow reusing the created scheme from models implemented in other programming languages. The scheme can be used on either CPU or GPU, and is distributed under the GPL3 license. Here, the LMC scheme is tested in a simple 0-dimensional adiabatic parcel model and then used in a 2-dimensional prescribed flow framework. The results are discussed with the focus on changes to the CCN sizes and compared 10 with other model simulations discussed in the literature. Copyright statement. This work is distributed under the Creative Commons Attribution 4.0 License


Introduction
libcloudph++ is an open-source library of schemes for representing cloud microphysics in numerical models.It was first introduced in Arabas et al. (2015) where the authors present the different microphysics schemes available in the library, show its user interface, and discuss its performance.The flagship component of libcloudph++ is the Lagrangian (i.e.particle tracking or particle based) microphysics scheme.The scheme resolves the evolution of aerosol, cloud droplet, and rain drop 1 size spectrum.
It allows representing from the first principles cloud microphysical processes and is especially well suited to track changes in the CCN size distribution that are caused by clouds (i.e.cloud-aerosol processing).The scheme can be used in models of any dimensionality or dynamical core, and can be run on both CPU and GPU.The open-source availability of the code, its clearly defined user interface, and the separation of concerns employed when designing libcloudph++ core code enable its usage and further development.
cloud-processed aerosol particles may again serve as CCN and influence microphysical properties of the next generation of clouds.The study of Pruppacher and Jaenicke (1995) estimates that on global average an atmospheric aerosol particle has been cycled 3 times by cloud systems.The LMC scheme introduced here offers a chance to represent the effects of such cloudprocessing on CCN sizes stemming from both collisions between water drops and aqueous phase oxidation reaction within water drops.To the authors knowledge, the presented scheme is the first to represent the impact of both collisions and aqueous phase chemistry on the aerosol size spectrum in the Lagrangian microphysics framework.
The structure of the presented work is as follows: Section 2 presents briefly the Lagrangian scheme available in lib-cloudph++.Section 3 discusses the design and user interface of the new aqueous chemistry scheme.Section 4 presents the validation tests of the new scheme.Section 5 discusses the results from simulations where the LMC scheme is incorporated into a simple model of a stratocumulus cloud.The effects of both collisions between water drops and aqueous phase oxidation of sulfur on the aerosol particle size distribution are presented.

Lagrangian microphysics scheme
The Lagrangian scheme used in this work is described in detail in Arabas et al. (2015) and this section only briefly summarizes its major concepts.In the Lagrangian approach to modeling cloud microphysics the computational domain is filled with "numerical point particles" representing a multiplicity of real particles (aerosol particles, cloud droplets or rain drops) of the same properties.Following the nomenclature introduced by Shima et al. (2009), the "numerical particles" are labeled here as superdroplets (SD).Each SD has a set of attributes describing the properties of the aerosol particles or water drops it represents.
As discussed in Arabas et al. (2015), for microphysical purposes, the required attributes are: multiplicity (N ), the position of SD in the computational domain, the wet radius (r w ), the dry radius (r d )2 and the hygroscopicity parameter (κ) 3 .The aqueous phase chemistry scheme extends the list of required attributes by masses of chemical compounds dissolved in droplets, see Sec. 3.
The key attribute of the Lagrangian microphysics scheme is the SD multiplicity.Multiplicity defines the number of aerosol particles or water drops represented by a given SD.All particles represented by one SD are assumed to be identical.Using multiplicites allows to reduce the complexity of the problem and enables efficient numerical computations.
The Lagrangian scheme used here requires no division into artificial categories of aerosol particles, cloud and rain water, as it is often done in the bulk schemes, for example Kessler (1995), Seifert and Beheng (2001), Morrison and Grabowski (2007).All the modeled microphysical processes are represented by calculating the changes to the SD attributes.The equation of condensational growth is solved iteratively for each SDs wet radius (see Sec. 5.1.3 in Arabas et al. (2015) for details).
The process of condensational growth from deliquescent aerosol particles to cloud droplets is thus resolved and no additional parametrisation of cloud droplet activation is required, as it is again often done in the bulk microphysics schemes (for example Morrison and Grabowski (2007)).Following Shima et al. (2009) the collisions between SDs are represented using Monte-Carlo scheme (see Sec. 5.1.4in Arabas et al. (2015) for details).
Information about SD attributes is retained within the model throughout the whole simulation.This means that the size distribution of both water drops and aerosol particles in each computational grid-cell can be easily obtained by taking into account the SD attributes of r w , r d and N .As a result, the Lagrangian scheme is capable of resolving the changes to both aerosol and water drop size distributions.The same functionality is offered by the 2-dimensional bin schemes, for example Ovchinnikov and Easter (2010) or Lebo and Seinfeld (2011).However, the Lagrangian approach greatly reduces the numerical diffusion errors.As discussed in Unterstrasser et al. (2016), it does introduce statistical errors, i.e. fluctuations between different realizations of the same collision/coalescence scenario.These errors are easier to minimize than diffusion numerical errors, for example by increasing the number of SDs in the computational domain or by averaging over an ensemble of simulation runs.
Lagrangian schemes solve ordinary differential equations instead of the partial differential equations encountered in the bin schemes, which is computationally more efficient.

Aqueous phase chemistry scheme
In order to represent the chemical composition of water drops the aqueous-phase chemistry scheme extends the list of SD attributes.The additional attributes are defined as the total mass of each of the chemical compounds in a given SD (including both the dissolved and, if applicable, dissociated fraction).An additional variable -the mass of the H + ions -is also added, in order to keep track of SD's pH.This results in eight new SD attributes needed for simulations with aqueous phase chemistry: the total mass of dissolved CO 2 (including: CO 2 * H 2 O, HCO 3 − and CO 2− 3 ), the total mass of dissolved NH 3 (including: the total mass of dissolved HNO 3 (including: HNO 3 (aq) and NO − 3 ), the total mass of created H 2 SO 4 (including: HSO − 4 and SO 2− 4 ) the total mass of H + ions.
The scheme needs to be coupled to the driver model that provides information about the environment in which SD are immersed (i.e.temperature, humidity, trace gas mixing ratio, and wind speed).The representation of aqueous phase chemistry more than doubles the number of required SD attributes and significantly increases the computational time.On the other hand, thanks to the added attributes, the mass of any ion for any SD can be easily diagnosed using just a dissociation constant.This, in turn, allows for a very straightforward representation of the aqueous chemical processes and does not call for any additional parametrisation.
All aqueous-phase chemistry included in the scheme is formulated under the assumption that solution droplets are diluted.
Therefore, in the LMC scheme, chemical processes are only performed for SDs with ionic strength smaller than 0.02 moles/liter (the same criterion is used for example in Ovchinnikov and Easter, 2010).In practice, this condition results in excluding from aqueous chemistry calculations SDs with small wet radii (i.e.SDs representing haze particles and very small cloud droplets).
Exclusion of SDs with small wet radii also prevents numerical issues during condensation procedure when changes in dry radius caused by oxidation could prevent convergence of the condensation scheme during the initial rapid growth of cloud droplets during activation.
Combining the Lagrangian microphysics scheme with aqueous phase chemistry is straightforward.Condensation/evaporation does not affect the chemical attributes of SDs.During collisions the mass of chemical compounds is summed when recalculating SD attributes (it is an extensive parameter).In principle the κ attribute should be recalculated in every time-step based on the new chemical composition of each SD.However, the κ values relevant for this study are very similar -the κ value of ammonium bisulfate is 0.61 (Petters and Kreidenweis, 2007) and of sulfuric acid is 0.64 (Kim et al., 2016).Therefore, the hygroscopicity parameter is assumed to be constant.
As discussed in Kreidenweis et al. (2003), one of the major sources of errors in modeling the in-cloud aqueous phase oxidation are the uncertainties when resolving the cloud droplet size distribution.It follows that combining the aqueous phase chemistry model with the very detailed Lagrangian microphysics scheme is beneficial.This will be discussed further in Sec. 4.

Dissociation
Dissociation is a reversible process of splitting of the molecules dissolved in water drops into ions.It is treated as an equilibrium process and is described using the dissociation constants.The dissociation constant of chemical compound A is denoted here by A .The dissociation constants are corrected for temperature using the formula of van 't Hoff (1885): where ∆H D denotes the reaction enthalpy of dissociation at constant temperature and pressure.The list of considered dissociation constants as well as their temperature correction coefficients is available in Tab.C2.The dissociation of water, although very small, is also taken into account 4 .No temperature correction is applied to water dissociation constant. 4The concentration of undissociated water molecules is so big that it is usually assumed constant and it traditionally multiplies the dissociation constant . (2) The [ ] brackets denote the concentration of each of the chemical compounds, (traditionally defined in units of moles per liter), capital letters denote chemical compound and roman numbers mark its oxidation state.In Eq. ( 2) the dissociation constants of SO 2− 3 , CO 2− 3 , and SO 2− 4 ions (i.e.HSO3 , HCO3 , and H2SO4 ) are multiplied by a factor of two, to take into account bigger electric charge number of those ions.Equation .(2)has only one unknown variable -the new equilibrium concentration of the H + ions.The new concentration is obtained iteratively using bisection algorithm 5 .The algorithm searches for solution between pH = -1 and pH = 9.The lower bound for the pH scale is unrealistically acidic and is only necessary at the start of the simulation when the initial SDs have very small volume and thus highly acidic pH.The upper bound is set arbitrarily, but is sufficient for the expected pH of the modeled droplets.At the end of the dissociation procedure the mass of H + ions is updated based on the new equilibrium concentration.
When the SD wet radius is quickly changing, for example during the initial condensational growth of cloud droplet or rain drop evaporation, the dissociation procedure requires small time-steps to reach convergence.The time-step used in dissociation procedure can be divided into user-specified number of sub-steps in order to prevent limiting the overall simulation time-step by dissociation.

Dissolution
The amount of the chemical compound that can dissolve into water drop from the gas phase is proportional to its partial pressure above the surface of the drop.Due to the longer timescale of the process, in contrast to dissociation, the transfer between the gas and liquid phase is not treated as an instantaneous process.Assuming that the water drop is internally mixed, the gas-liquid transfer is limited by the diffusion of gas phase particles to the drop surface (gas-phase limitation) and the probability that the molecule will enter the drop after collision (interfacial limitation).Following chapter 8.4.2 in Warneck (1999), for a chemical compound "A" the rate of transfer from the gas phase to the aqueous phase equals is the effective equilibrium constant for dissolution of the chemical compound "A".The dissolution constants depend on the temperature following similar relation as for dissociation Eq. ( 1).Table C3 shows the equilibrium dissolution constants and their temperature corrections and Tab.C4 presents the diffusion and mass accommodation coefficients.The term "effective" marks that the dissolution constants take into account the increase of the efficiency due to dissociation (see Seinfeld and Pandis (1998) chap. 7.4 for the exact equations).Equation ( 3) is solved for each SD and for each of the considered trace gases.It is solved implicitly with respect to aqueous-phase concentration and explicitly with respect to the gas-phase concentration.The input ambient trace gas concentration is calculated from the trace gas mixing ratio provided by the driver model to which the LMC scheme is coupled.
Obtained aqueous phase concentration is recalculated to the mass of dissolved chemical compounds and the corresponding SD attribute is updated.The changes in the ambient trace gas mixing ratios are calculated by LMC scheme by summing the changes in chemical composition in all SDs in a given grid-cell and then subtracting them from the trace gas mixing ratio of the driver model.This approach does not ensure that per each time-step the total dissolved mass of a given trace gas does not exceed the available ambient mixing ratio.To prevent that, relatively short time-steps should be applied.If necessary the user can divide the model time-step into sub-steps.

Oxidation
The reaction rates of oxidation by ozone and hydrogen peroxide can be described as (Hoffmann and Calvert, 1985): where A is the reaction rate of the chemical compound "A" and k 0,...,4 are the reaction rate coefficients.k 0,...,4 depend on the temperature following similar relation as for dissociation Eq. (1).Table C5 shows the values of reaction rate coefficients and their temperature correction coefficients.
Equations ( 4) and (5) return the new concentration of H 2 SO 4 created in each SD in each time-setp.Based on the new concentration, the new mass of H 2 SO 4 and the new dry radius are calculated and the corresponding SD attributes are updated.
The dry particle density of 1.8 g/cm 3 is assumed while evaluating the dry radius from the H 2 SO 4 mass.
For the typical atmospheric conditions, say pH between 3 and 6 (i.e.

Initialization
The initial aerosol is assumed to be ammonium bisulfate (NH 4 HSO 4 ), with dry particle density of 1.8 g cm −3 .Using dry particle density and dry radius of each SD, the initial mass of H + , NH + and SO 2− 4 ions is calculated.The initial mass of other molecules and ions is equal zero and is therefore not in equilibrium with the initial ambient trace gas conditions.For the initial conditions above supersaturation it is advisable to allow for a spin-up period with only condensation/evaporation and the equilibrium chemical processes enabled, to allow the model to reach equilibrium.

Validation
The LMC scheme is set to reproduce results from model intercomparison study by Kreidenweis et al. (2003), where several bulk and bin schemes representing cloud microphysics and aqueous-phase chemistry were tested in an adiabatic parcel model setup.A parcel model is a 0-dimensional model that represents an idealized scenario of a finite volume of air rising adiabatically with a constant vertical velocity.As the parcel of air raises, its temperature decreases leading to supersaturation.This results in activation and further condensational growth of cloud droplets.For the studied oxidation reaction, the presence of liquid water enables aqueous-phase chemical reactions and leads to creation of sulfuric acid within cloud droplets.The collisions between cloud droplets are not included in the parcel simulations to allow an easy comparison with Kreidenweis et al. (2003).
The initial conditions are the same as in Kreidenweis et al. (2003) and are provided for convenience in Tab. 1.The simulation starts below cloud base (i.e. with subsaturation).The initial aerosol is ammonium bisulfate and the initial aerosol particle size distribution is assumed to be lognormal with one mode where n(r d ) is the spectral density function of aerosol particle sizes, n tot is the total aerosol concentration, r d is median radius and σ g is the geometric standard deviation.
The parcel model employed in this study uses dry air density ρ d , dry air potential temperature θ, water vapor mixing ratio r v , and mixing ratios of ambient trace gases as model variables.In order to calculate ρ d at each time-level (or each height-level of the parcel ascent) the model assumes a vertical profile of pressure.In the presented simulations the assumed pressure profile is obtained by integrating the hydrostatic equation and assuming that density of air is constant at each height-level (piecewise constant).Then, at each level, ρ d is calculated from the ideal gas law taking into account the current r v and θ.Because the simulated air parcel is assumed to be adiabatic, only processes resolved by the Lagrangian scheme can change θ, r v and other trace gas mixing ratios.In each model time-step, the Lagrangian microphysics scheme changes θ and r v according to Eq. 25 and 26 from Arabas et al. (2015).The changes in the trace gas mixing ratios are resolved following procedure discussed in Sec.3.2.It is assumed that the initial mass of dry air within the parcel is 1 kg.
Figure 1 shows the general physical and chemical conditions from the cloud base up to the end of the test run 1.2 km above the cloud base.Two vertical axes are used, representing either the time or the height above the cloud base.Figure 1a shows the liquid water mixing ratio (LWC).The increase in LWC is linear and the LWC reaches above 2 g/kg at height 1.2km above the cloud base.Figure 1b shows the total SO 2 concentration (both in gas phase and dissolved in water) in ppb units.The concentration of SO 2 is decreasing due to oxidation taking place in cloud droplets.Figure 1c shows the water weighted average pH of cloud droplets.The pH near the cloud base is very low due to acidic nature of the assumed initial aerosol and small size of activated cloud droplets.As the drops grow bigger and become more diluted, the average pH increases.Figure 1 compares well with Fig. 1 from Kreidenweis et al. (2003), the overall differences between the two figures are small.The biggest discrepancy 5 is in the LWC profile.This is arguably caused by different pressure profiles in the two parcel models (the pressure profile in Kreidenweis et al. (2003) is not specified).
At the end of the test simulation, 84% of SO 2 is converted into H 2 SO 4 and the final water weighted average pH is equal to by bin schemes.
Microphysics schemes taking part in the Kreidenweis et al. (2003) intercomparison study reported significant differences between the number of activated cloud droplets.Based on Tab. 1 in Kreidenweis et al. (2003) the droplet number concentration at the cloud base varied between 275 and 358 cm −3 .One of the differences between the bin schemes responsible for causing this discrepancy is the different water vapor mass accommodation coefficient α M leading to different predicted maximum supersaturation.Figure 6 in Kreidenweis et al. (2003) shows that the observed maximum supersaturations were lower (between 0.23% -0.26%) for schemes using high values of α M (either 0.5 or 1).In contrast, a scheme using α M = 0.042 predicted maximum supersaturation equal to 0.37%.The Lagrangian scheme used in this study reports concentration 264 cm −3 at RH = 1 and 272 cm −3 at the level of maximum supersaturation.The maximum supersaturation is equal to 0.27%.The Lagrangian scheme assumes α M equal to unity and therefore fits with the trend of high α M causing lower supersaturation presented in Fig. 6 in Kreidenweis et al. (2003).
Another cause for the discrepancy between the bin schemes listed in the intercomparison study are the different sizes and locations of bins in different models.Along those lines, here it is tested how sensitive the Lagrangian scheme is to the number of SDs.The results of this test are summarized in Fig. 2 showing the cloud droplet concentration at the cloud base (a), the maximum supersaturation (b), the average pH (c) and the total sulfate production (d).The results are plotted against the logarithm of base two of the number of SDs in the computational domain (meaning that "0" represents one SD and "10" little (between 265 cm −3 and 274 cm −3 ).The cloud droplet concentration from simulation with 128 SDs is the largest outlier.
The concentrations from simulations with SD number between 256 and 1024 vary between 274 cm −3 and 272 cm −3 .For SD numbers between 32 and 128 there are insignificant changes in the maximum supersaturation.The values of pH vary by 0.01 and the total sulfate production increases by 1 ppt.There are, however, large differences between the number of droplets at the cloud base (between 265 cm −3 and 335 cm −3 ).This confirms the observations from Kreidenweis et al. (2003) that the predicted  The effect of in-cloud sulfate production on aerosol particle size distribution presented in Fig. 3, combined with other tests presented in this chapter, document the correctness of the implementation of aqueous chemistry in the Lagrangian scheme.
The formation of the "Hoppel minimum" was reported by many observational studies, see Hoppel et al. (1994); Bower et al. (1997); Hudson et al. (2015).Figure 3 compares well with aerosol size distribution plots from the intercomparison study shown in Fig. 9 in Kreidenweis et al. (2003).Other numerical schemes also reported the formation of "Hoppel minimum", see for example Flossmann (1994); Feingold andKreidenweis (2000, 2002); Ovchinnikov and Easter (2010).The collisions between water drops are represented using the geometric kernel with collision efficiency for big drops from Hall (1980) and for small droplets from Pinsky et al. (2008).For big drops, the collision efficiencies were obtained from the fit to measurements, see Hall (1980).For small droplets, the collision efficiencies were based on numerical simulations taking into account turbulence typical for stratocumulus clouds, see Pinsky et al. (2008).The collision efficiencies are provided via a look-up tables for different drop sizes.
The initial conditions are summarized in Tab. 2. The computational domain size is 1.5 km in both directions and the computational grid is composed of 75×75 cells of equal size (the grid lengths are 20 m) and is periodic in the horizontal direction.
The initial air density profile corresponds to the hydrostatic equilibrium with the pressure of 1015 hPa at the bottom of the domain.At the beginning of the simulation it is assumed that there is no condensed water, and the initial profiles of θ and r v are constant with altitude.To keep the simulation setup simple and due to a relatively low vertical extent of the computational domain, the initial trace gas volume fractions are also assumed to be constant with altitude.This unrealistic initial condition results in very high initial supersaturation in the upper part of the domain.As a consequence a 10 5 second (∼ 2h 45min) spin-up period is necessary to allow for the simulated water drops to reach equilibrium with their environment.During the spin-up only the reversible processes (condensation and evaporation, dissolving of trace gases and dissociation into ions) are allowed and the supersaturation is limited to 5% (RH=1.05).After spin-up the simulations are run for 30 minutes.The chosen simulation time is enough to deplete the SO 2 available in the cloudy part of the domain as well as to create precipitation.
Similarly to the adiabatic parcel test, the initial aerosol is ammonium bisulfate and the aerosol particle size distribution is lognormal with one mode.The initial condition for trace gases is defined in terms of volume fractions and then translated to mixing ratios that serve as the the model variables.The initial SO 2 , O 3 and H 2 O 2 volume fractions are taken from the simulation setup used in Ovchinnikov and Easter (2010).The values for SO 2 and O 3 are based on the measurements from MASE campaign (Wang et al., 2008) and the value for H 2 O 2 is based on the representative values for the Eastern Pacific Ocean (Genfa et al., 1999).The NH 3 , HNO 3 and CO 2 volume fractions are the same as in the parcel test from Sec. 4.
The setup detailed in Tab. 2 corresponds to "very clean conditions" (i.e.low aerosol particle concentrations).Initial aerosol particle sizes are also relatively small.Three additional simulation cases are studied to check the sensitivity of the model to different conditions.In case1 the reversible chemical processes are allowed, but oxidation is prohibited.In case2 the initial volume fraction of NH 3 is increased and in case3 the initial aerosol size distribution is changed.The conditions for all the sensitivity simulation cases are summarized in Tab. 3.
As discussed in Flossmann (1994), the initial chemical scenario is idealized.For instance, although the initial conditions represent clean maritime environment, the setup lacks sea salt aerosol particles.As discussed by Twohy et al. (1989)  to include aqueous chemistry into the Lagrangian cloud microphysics and therefore the decision was made to start with the simplified setup.
The initial aerosol size distribution parameters are based on the test cases studied in Feingold and Kreidenweis (2002).
The discussion presented in their study introduced two regimes for oxidation with regard to the mean aerosol size r d and precipitation: (i) for relatively small initial r d production of sulfate enhances precipitation, (ii) for relatively big initial r d production of sulfate suppresses precipitation.The overall impact depends strongly on the initial concentration of aerosol particles, see Feingold and Kreidenweis (2002) for the discussion.The short simulation time used in this study hinders analysis of the impact of oxidation on the overall precipitation.The work presented here focuses on the evolution of aerosol particle sizes and pH values of cloud and drizzle droplets.Future LES simulations should focus on the impacts of aqueous chemistry on precipitation, cloud lifetime, and cloud dynamics.
The kinematic setup precludes any links between cloud microphysical processes and dynamics of the air motion.The setup limits the study to the smooth velocity and therefore smooth saturation fields and prevents mixing between air parcels with different trajectories and properties.On the other hand, the kinematic setup has low computational cost and allows easy testing and sensitivity analysis.Prescribing the velocity ensures that all changes to the aerosol particle and water drop size distributions are caused by the cloud microphysics and aqueous-phase chemistry alone.Moreover, the kinematic setup allows for a straightforward selection of the updraft and downdraft regions, further simplifying the analysis of the microphysical processes.

Results
Figure 4 shows the model state after 30 minutes of simulation from base case (see Tab. 2). Figure 4a shows the concentration of unactivated aerosol particles (defined as the SDs with wet radius smaller than 1 µm).The lower part of the plot (below 900 m) shows cloud-free conditions and corresponds to the initial concentration of aerosol particles.The upper part of the plot shows the interstitial aerosol particles, i.e. those aerosol particles that did not activate.The difference between the upper and lower parts of Fig. 4a shows the impact of nucleation scavenging on aerosol population.The regions with slightly higher concentration of the in-cloud aerosol particles near the cloud base correspond to regions with low vertical velocities, lower supersaturations and thus lower concentrations of cloud droplets.Figure 4b shows the concentration of cloud droplets (defined as the SDs with wet radii between 1 and 25 µm).It is nearly constant with height, that agrees with the observations in stratocumulus clouds (e.g.Pawlowska et al., 2000).The regions with lower cloud droplet concentrations correspond to the regions with drizzle (see Fig. 4f). Figure 4c shows the rain water mixing ratio (water drops with wet radius greater than 25 µm) using a logarithmic color scale.Rain forms quickly in the simulation due to the relatively high values of cloud droplet radii after the spin-up caused by the low initial aerosol particle concentration.The footprint of precipitation can be seen in Fig. 4a and f where the cloud droplet concentration is depleted in regions of drizzle.Figure 4d shows the mean dry radius of all particles (both aerosol particles and water drops).The mean dry radius is increasing due to oxidation.and concentration of S VI molecules (f).The thresholds for particle radii are: aerosol < 1 µm; 1 µm < cloud < 25 µm; rain > 25 µm.Note the logarithmic scale for the rain water mixing ratio plot.
particle size distribution.Figure 4e depicts the cloud droplet effective radius.As expected, the effective radius increases with height.At the top of the cloud the effective radius reaches 20 µm, which is linked to the small cloud droplet concentration.
High effective radii imply efficient drizzle production after the spin-up (usually water drop radius ∼ 12 µm is reported as the threshold value for efficient collisions between water drops and the production of precipitation, for example Rosenfeld and Gutman (1994); Pawlowska and Brenguier (2003)).1.The increase of pH with height is also observed in 1-dimensional model representing processing of sulfur in small cumuli in marine environment (Alfonso and Raga, 2002).Due to the pH variability shown in The study by Walcek and Taylor (1986) also reported that the pH of droplets increased with height due to dilution despite the production of sulfuric acid.In turn, increased pH promotes oxidation by O 3 in the upper parts of the cloud, whereas oxidation by H 2 O 2 dominates in lower parts of the cloud, according to their study.
Case1 shown in Fig. 5b represents a hypothetical "no oxidation" scenario where all physical and chemical conditions are the same as in base case, the reversible chemical processes are allowed, and oxidation is prohibited.The scenario without oxidation is overall less acidic than base case (Fig. 5a).Additionally, without oxidation there is no difference between the pH values in the updraft and downdraft in Fig. 5b.Without oxidation, all the chemical processes are reversible and the dissolved chemical compounds are outgassed to the atmosphere as the cloud droplets evaporate in the downdraft.
Case2 differs from base case by increasing the initial NH 3 volume fraction from 0.1 ppb-v to 0.4 ppb-v (see Tab. 2 and Tab.3).Because the initial aerosol particle size distribution is the same as in base case, the mean aerosol and droplet sizes and concentrations at the end of the simulation are not different from base case (not shown).Figure 5c shows the liquid water weighted average pH for case2.The average pH in case2 (Fig. 5c) is higher than in base case (Fig. 5a), that is, both cloud droplets and rain drops are less acidic in case2 than in base case.In contrast to base case, in the updraft (left-hand side of the plots), the pH in case2 actually decreases with height above the cloud base.This is because the higher initial NH 3 volume fraction increases its uptake and counters the low pH values caused by initial acidic aerosol particles.Then, as the water drops are advected upwards, oxidation produces sulfuric acid and the average pH decreases.Near the cloud top, the NH 3 is degassed back to the environment.Case2 results are in agreement with the trajectory ensemble model simulations by Zhang et al. (1999).In their study, the initial aerosol size distribution is the same as in base case and case2.However, their initial trace gas volume fractions are much higher and aim to represent a "moderately polluted marine environment" (their base case NH 3 volume fraction is ten times larger than base case value assumed here).As in case2 presented here, the high initial NH 3 volume fractions in Zhang et al. (1999) increase the pH near the cloud base and promote oxidation by O 3 during the first minutes after the simulated parcels entered the cloud.Because the sulfuric acid was produced, the pH dropped and oxidation by H 2 O 2 becomes dominant in the higher regions of the cloud, as reported in their study.
Case3 increases the initial aerosol concentration to 150 cm −3 , while keeping all other initial conditions the same as in base case (see Tab. 2 and Tab.3).In general, higher initial aerosol particle concentration results in higher cloud droplet concentrations.This in turn creates smaller cloud droplet effective radii that virtually prohibits the onset of precipitation during the 30 minutes simulation time (not shown).Figure 5d shows the liquid water weighted average pH for case3.Similar to base case (Fig. 5a), the pH increases with height due to the dilution and the downdraft droplets are more acidic due to the ongoing oxidation.However, case3 is more acidic than base case because the overall droplet sizes are smaller and they are therefore 18 At the end of base case simulation, 18% of the total available S IV is oxidized.As a result, 0.14 µg/m 3 of dry particulate matter are created during oxidation (an average value for the whole computational domain reported in relation to the dry air volume).In total, 40% of the final dry particulate matter is created due to oxidation and 60% originates from the initial aerosol mass.The oxidation is a significant source of dry particulate matter because the initial aerosol mass is very low (only 0.21 µg/m 3 dry air).Oxidation by H 2 O 2 is the dominant path: 92% of the S VI molecules created during oxidation are oxidized by H 2 O 2 .More alkaline conditions of case2 enhance the efficiency of oxidation.At the end of case2, 21% of available S IV is oxidized.As a result, oxidation produces 0.16 µg of dry particulate matter per m 3 of dry air (average over the whole computational domain).For case2, 44% of the final dry particulate matter is created due to oxidation and 56% originates from the initial ammonium bisulfate aerosol.Similarly to base case, the significance of oxidation as a source of dry particulate matter is caused by a very low initial aerosol mass.Due to more alkaline conditions, oxidation by O 3 becomes more important than in base case.At the end of case3 simulation, 39% of the S VI molecules that are created during oxidation are produced by the O 3 path and 61% by the H 2 O 2 path.In contrast, more acidic conditions of case3 hinder the O 3 reaction path.Virtually all molecules of sulfate that are created during oxidation are oxidized by H 2 O 2 .As a result, the conversion of sulfur to sulfate is slightly less effective in case3.At the end of case3 simulation, 17% of available S IV is oxidized.As a result 0.13 µg of dry particulate matter are created per m 3 of dry air.At the end of case3 simulation, 17% of the dry particulate matter is created by oxidation and 83% originates from the initial aerosol.The initial aerosol mass is larger in case3 than in base case due to the higher initial aerosol concentration (case3 contains initially 0.61 µg/m 3 of dry particulate matter).Because of this, even though the produced sulfate mass is only slightly lower than in base case, the relative importance of oxidation decreases by more than 20 percentage points in case3.Due to the simple kinematic setup chosen in this study the values reported here cannot be treated as representative for the atmospheric conditions.They are shown to allow comparison between base case and the sensitivity test cases.
Finally, the impact of collisions and aqueous phase oxidation of sulfur on the aerosol and water drop size distributions is examined.For this purpose, the aerosol particle size distributions from base case (Fig. 6a) and case3 (Fig. 6b) are compared.The black line represents the initial aerosol size distribution and the green and red lines represent the final aerosol size distribution for in-cloud (r c > 0.01 g/kg) and precipitating (r r > 0.01 g/kg) grid-cells.The two cases are chosen because they have different initial aerosol size distributions.In both cases the cloud-processed aerosol size distributions (green and red lines) have a bimodal shape.This is a footprint of oxidation that creates the Hoppel minimum in the dry radius size distribution.The same effect is obtained in the adiabatic parcel tests discussed in Sec. 4 and shown in Fig. 3.Moreover, the efficient collisions between water drops in base case create a tail of bigger aerosol sizes in Fig. 6a.The effect is stronger for the precipitating grid-cells (red line).In case3 fewer collisions between water drops occur than in base case and therefore no precipitation and no tail of big aerosol particles is created.Also, in case3, the change in size distribution of aerosol particles caused by oxidation is smaller because the produced sulfate is divided among bigger number of aerosol particles.

Summary and outlook
The work presented here describes a new extension of the libcloudph++ that allows including aqueous phase chemical reactions within water drops in the Lagrangian microphysics scheme.The extension covers the aqueous phase oxidation of sulfur to sulfate.The modular way in which the library is implemented along with the provided documentation should allow, if needed, further development to cover more chemical compounds and reactions.The 0-dimensional and 2-dimensional tests described in this work as well as comparison with other numerical studies using bin microphysics schemes along with aqueous chemistry representation document the correctness of the design and the implementation of the LMC scheme.Additionally, the changes in the user interface due to aqueous chemistry extension are described in Sec.A in the Appendix.Section C in the Appendix completes the description with a list of chemical constants used in the library and chemical reactions included.
The models used in this study to test the chemistry scheme provide a simplified view of the macrophysical cloud properties.
They enable validation and testing of the Lagrangian scheme but do not provide a good balance between the representation of cloud microphysics and dynamics.As a next step, the Lagrangian scheme needs to be coupled to an eddy-resolving model.
This would allow quantifying how microphysical and chemical processes affect precipitation in the model and how they affect the cloud lifetimes simulated by the model.
ssumed to be ammonium bisulfate particle density of 1.8 g cm 3 .sity and dry radius of each SD, the + and SO 4 2ions is calculated.ther molecules and ions is equal ot in equilibrium with the initial ditions.For the initial conditions it is advisable to allow for a spinondensation/evaporation and the processes enabled, to allow the rium.
f the Lagrangian microphysics ++ is presented in Sec 5.2. in Here, additional information reus phase chemistry scheme is pro-h++ is implemented in C++ and nclature related to this program-.For a thorough introduction to guage see Stroustrup (2013).istry module is implemented as to the Lagrangian microphysics to true at the end of spin-up.Other parameters that cannot be changed during simulation are encapsulated in lgrngn::opts_init_t structure.Chemistry module adds three fields to this structure: (i) A Boolean chem_switch field that enables memory allocation for additional variables needed for chemistry representation.(ii) An integer sstp_chem field that defines the number of sub-steps to be carried out in aqueous chemistry calculations.(iii) A real_t chem_rho field that defines the dry aerosol density, see Code Listing.3.2.The user interface of the Lagrangian microphysics scheme of libcloudph++ is presented in Sec 5.2. in Arabas et al. (2015).
Here, additional information related to the new aqueous phase chemistry scheme is provided.The libcloudph++ is implemented in C++ and therefore some nomenclature related to this programming language is used.For a thorough introduction to C++ programming language see Stroustrup (2013).
The aqueous chemistry module is implemented as an optional extension to the Lagrangian microphysics scheme in lib-10 cloudph++.It uses the same libcloudphxx::lgrngn namespace as the original scheme.Again the template parameter real_t selects between floating point formats of simulations.The Lagrangian microphysics scheme options are grouped into a structure named lgrngn::opts_t.Chemistry module adds three Boolean fields to this structure: chem_dsl chem_dsc and chem_rct, see code listing in Fig. A1.When set to true by the user, they switch on dissolving of trace gases into water drops, dissociation of chemical compounds in water drops and oxidation reaction, respectively.The parameters in lgrngn::opts_t can be changed 15 during simulation.For example during the 2-dimensional kinematic simulations from Sec. 5, oxidation is enabled by setting the chem_rct parameter to true at the end of spin-up.Other parameters that cannot be changed during simulation are encapsulated in lgrngn::opts_init_t structure.Chemistry module adds three fields to this structure: (i) A Boolean chem_switch field that enables memory allocation for additional variables needed for chemistry representation.(ii) An integer sstp_chem field that defines the number of sub-steps to be carried out in aqueous chemistry calculations.(iii) A real_t chem_rho field that defines  The init() method performs initialization and should be called first.As discussed in Arabas et al. (2015) the first three arguments are obligatory and should point to the dry air potential temperature, water vapor mixing ratio and dry air density fields of the driver model that uses the libcloudph++.The next three arguments should point to the Courant number field components.They are optional and depend on the dimensionality of the solved problem.For example, for the parcel model tests from Sec. ?? none are necessary, whereas for the 2-dimensional kinematic model from Sec. 5 two arguments are specified in order to describe the velocity field.The last argument of init() is a map with keys from chem_species_t enumerator and values pointing to the corresponding trace gas mixing ratio fields from the driver model.This is an optional argument for The init() method performs initialization and should be called first.As discussed in Arabas et al. (2015) the first three arguments are obligatory and should point to the dry air potential temperature, water vapor mixing ratio and dry air density fields of the driver model that uses the libcloudph++.The next three arguments should point to the Courant number field components.They are optional and depend on the dimensionality of the solved problem.For example, for the parcel model tests from Sec. ?? none are necessary, whereas for the 2-dimensional kinematic model from Sec. 5 two arguments are specified in order to describe the velocity field.The last argument of init() is a map with keys from chem_species_t enumerator and values pointing to the corresponding trace gas mixing ratio fields from the driver model.This is an optional argument for simulations with aqueous phase chemistry.
During time-stepping the Lagrangian scheme computations are performed by step_sync() and step_async() methods.The first one gathers all the processes that a ect the driver model fields (such as condensation/evaporation or aqueous phase chemistry) and the second one gathers all th be calculated asynchronously (for sedimentation).The list of argum method is extended by the chemis to the init() method, a map link enumerator items with the driver fields needs to be provided as the la The Lagrangian scheme overwrit fields during simulation.The sig method is not changed by the new As discussed in Arabas e lgrngn::particles_t structure pro for obtaining statistical information parameters (prefixed with diag).adds to them the diag_chem met  and the second one gathers all the processes that can be calculated asynchronously (for example collisions or sedimentation).The list of arguments of step_sync() method is extended by the chemistry module.Similar to the init() method, a map linking chem_species_t enumerator items with the driver model mixing ratio istry) and the second one gathers all the processes that can be calculated asynchronously (for example collisions or sedimentation).The list of arguments of step_sync() method is extended by the chemistry module.Similar to the init() method, a map linking chem_species_t enumerator items with the driver model mixing ratio fields needs to be provided as the last optional argument.The Lagrangian scheme overwrites the driver model fields during simulation.The signature of step_async method is not changed by the new chemistry module.C2.Dissociation constants and their temperature correction coefficients (taken from Kreidenweis et al., 2003).Dissociation of H 2 SO 4 is taken from Tab. 6.A.1 in Seinfeld and Pandis (1998).Table C3.Dissolution constants and their temperature correction coefficients (taken from Kreidenweis et al., 2003).
oxidation reaction path reaction rate coefficient (liter moles -1 s -1 ) at 298K temperature correction −E R (K) of water.This leads to a different definition of the dissociation constant for water: H2O = [H + ][OH − ] 5 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-96Manuscript under review for journal Geosci.Model Dev. Discussion started: 23 April 2018 c Author(s) 2018.CC BY 4.0 License.It is assumed that there is no electric charge of water drops and therefore the concentrations of positive and negative ions created during dissociation should balance each other.Using dissociation constants (see Tab. C2), all ion concentrations can be expressed as a function of the total concentration of the dissolved chemical compounds and the concentration of H + ions.The neutral charge condition can be expressed as positive ions algorithm from Boost library.See www.boost.org/docfor documentation and Alefeld et al. (1995) for derivation 6 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-96Manuscript under review for journal Geosci.Model Dev. Discussion started: 23 April 2018 c Author(s) 2018.CC BY 4.0 License.where D A and α M A are the diffusion and mass accommodation coefficients for the chemical compound "A", <v> = 8RT π M A is the average velocity of the molecules calculated from the Maxwell-Boltzmann distribution function, M A is the molar mass of the chemical compound "A", c A∞ is the ambient concentration of the trace gas "A" and e f f A

4. 83 .Figure 1 .
Figure 1.Physical and chemical conditions in the adiabatic parcel model.Panel a shows the liquid water mixing ratio (LWC), panel b shows the SO 2 concentration (both gas phase and dissolved) panel c shows the water weighted average pH of the simulated population of water drops.

Figure 2 .
Figure 2. Results of the convergence test for the adiabatic parcel simulations.All figures show how a given parameter depends on the number of SDs (shown on the abscissa as the logarithm of base two of the number of SDs).Panel a the cloud droplet concentration at the cloud base, panel b the maximum supersaturation, panel c the water weighted average pH at the end of simulation and panel d the total sulfate production.

5
cloud droplet number concentration strongly depends on the representation of the size distribution of modeled aerosol particles and cloud droplets and that this may become a major source of uncertainties in the microphysics representation.Decrease in the SD number below 32 leads to a big variance in the cloud droplet concentration as well as other parameters.

Figure 3 Figure 3 .
Figure 3 shows the simulated modification of the aerosol size distribution.Red line depicts the initial distribution and the green line shows model state at the end of adiabatic parcel test.For convenience, Fig. 3 uses both logarithmic (left panel) and 10 kinematic modelThe kinematic model mimics a single 2D eddy spanning a stratocumulus cloud deck and a boundary layer below it.The model is based on a test scenario from the 8 th International Cloud Modeling Workshop (ICMW;Muhlbauer et al., 2013, case 1).The velocity field is prescribed as inSzumowski et al. (1998);Morrison and Grabowski (2007);Rasinski et al. (2011).The same model was used when presenting the initial release of libcloudph++, see Sec. 2 inArabas et al. (2015) for the details of the model formulation.The model operates on the Eulerian grid.At each model time-step, the temperature, moisture, and trace gas fields are advected using the prescribed velocity field.Then, the model variables are passed to the LMC scheme, where the microphysical and chemical processes are resolved.Finally, the source and sink terms due to microphysics and chemistry are 12 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-96Manuscript under review for journal Geosci.Model Dev. Discussion started: 23 April 2018 c Author(s) 2018.CC BY 4.0 License.calculated and applied in each model grid-cell as described in Sec. 2 and 3.
, sea salt aerosol particles are alkaline, which may in turn increase the pH of water drops and thus affect the oxidation rate.On the other hand, a study by vonGlasow and Sander (2001) indicates that alkaline sea salt particles are quickly converted to acidic due to the uptake of HCl vapor.More importantly, including sea salt would result in aerosol particles with very different hygroscopicity values (κ of ammonium bisulfate is 0.61, whereas κ of NaCl is 1.28;Petters and Kreidenweis, 2007).Including 13 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-96Manuscript under review for journal Geosci.Model Dev. Discussion started: 23 April 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 4 .
Figure 4. Base case setup (see Tab. 2).All panels depict model state after 30 minutes simulation time (excluding the spin-up) and show: aerosol concentration (a), cloud droplet concentration (b), rain water mixing ratio (c), mean dry radius (d), cloud droplet effective radius (e) Figure 4f shows the concentration of S VI molecules (all molecules containing sulfur at +6 oxidation state) and represents molecules from the initial ammonium bisulfate (NH 4 HSO 4 ) aerosol and the molecules created during oxidation.It corresponds to the mean dry radius plotted in Fig. 4d.Additionally, some effects of collisions and precipitation can be seen when comparing the irregular features from Fig. 4f with rain water mixing ratio in Fig. 4c.Precipitation displaces the largest water drops, which causes the irregular distribution of S VI molecules in cloudy 16 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-96Manuscript under review for journal Geosci.Model Dev. Discussion started: 23 April 2018 c Author(s) 2018.CC BY 4.0 License.grid-cells.Figure 4f also shows that the Lagrangian scheme can track the dissolved chemical compounds in the evaporating rain drops below the cloud base.

Figure5Figure 5 .
Figure 4b, d, e and f show a layer of very clean air above the cloud which is caused by sedimentation of cloud droplets.In the downdraft region, the prescribed velocity field advects the clean layer into the domain.This feature is not present in the aerosol concentration plot (Fig. 4a) because the clean layer contains small aerosol particles with small sedimentation velocity.5

Figure 5
Figure5shows the liquid water weighted average pH in each computational grid-cell from base case (a) and sensitivity test cases (b-d).In order to better adjust the color scale to the in-cloud pH variability, pH values below 3 that correspond to very 10

Fig
Fig. 5a oxidation by O 3 happens mostly near the cloud top in base case.As discussed in Sec.3.3, the rate of oxidation by O 3 increases significantly with increasing pH (see Eq. 4) whereas oxidation by H 2 O 2 does not depend on the acidity (Eq.5).

19Figure 6 .
Figure 6.Size distributions of dry radii for base case (a) and case3 (b).The initial dry radius size distribution is marked in black, final dry radius size distribution from grid-cells with r c > 0.01 g/kg in green and from grid-cells with r r > 0.01 g/kg in red.See Tab. 2 and Tab. 3 for a definition of simulation setups.

20
the dry aerosol density, see code listing in Fig. A2.The names of chemical compounds available in the aqueous phase chemistry module are stored in a chem_species_t enumerator, see code listing in Fig. A3.The state of all variables used by the Lagrangian scheme is stored in an instance of the lgrngn::particles_t structure shown in code listing in Fig. A4.The second template parameter of that structure selects between 25 21 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-96Manuscript under review for journal Geosci.Model Dev. Discussion started: 23 April 2018 c Author(s) 2018.CC BY 4.0 License.
[H +] between 10 −3 and 10 −6 ), it can be said that the rate of oxidation by H 2 O 2 does not depend on pH (see Tab. C2 for the dissociation constant values).In contrast, oxidation by ozone depends strongly on pH of the solution and can become very fast if pH is high.For example, increasing pH by 1 point results in approximately 100 increase in O 3 reaction rate.

Table 1 .
Initial conditions for the adiabatic parcel test.

Table 2 .
Initial conditions for base case of 2-dimensional kinematic model.

Table 3 .
Initial conditions for sensitivity test cases of 2-dimensional kinematic model.Specified are: aqueous phase chemistry choice, initial volume fraction of NH 3 , mean radius of the assumed lognormal aerosol particle size distribution r d , total aerosol concentration n tot and geometric standard deviation σ g .Other parameters for each case are the same as in base case (Tab.2) The parameters that distinguish each sensitivity test case are marked in bold.case oxidation reaction NH 3 [ppb-v] r d [µm] n tot [cm −3 ] σ g different condensational growth of aerosol particles.The setup used in this study also lacks other particles containing sulfate, such as ammonium sulfate or sulfuric acid aerosol particles.The reason behind the chosen setup, is that this is the first attempt 14 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-96Manuscript under review for journal Geosci.Model Dev. Discussion started: 23 April 2018 c Author(s) 2018.CC BY 4.0 License.

Table C1 .
Chemical compounds considered in this work.

Table
for relevant chemical compounds .diffusion coeff.D A (m 2 /s) mass accommodation coeff.α M A