libcloudph++ 2.0: aqueous-phase chemistry extension of the particle-based 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 particle-based microphysics scheme with a Monte Carlo coalescence 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 particle-based microphysics and chemistry scheme allows for tracking of the changes in the cloud condensation nuclei (CCN) distribution caused by both collisions between cloud droplets and aqueous-phase oxidation. The scheme is implemented in C++ and equipped with bindings to Python. The scheme can be used on either a CPU or a GPU, and is distributed under the GPLv3 license. Here, the particle-based microphysics and chemistry 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 a focus on changes to the CCN sizes and comparison with other model simulations discussed in the literature.


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 programming interface, and discuss its performance. The flagship component of libcloudph++ is the particle-based (i.e. parti- 15 cle tracking or "Lagrangian-in-droplet-radius-and-space") 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 main software design principle employed while developing libcloudph++ core code is the separation of concerns. The code is open- 20 source and it's programming interface is documented in Arabas et al. (2015). All those features facilitate further development 1 For convenience cloud droplets and rain drops will be often labeled together as water drops 1 and usage of libcloudph++.
Sulfate aerosols cool Earth's climate by scattering sunlight and thus increasing Earth's shortwave albedo (direct radiative forcing) and by changing radiative properties of clouds (cloud albedo effect). According to the chapter 8 of IPCC 2 Assessment Report (Myhre et al., 2013), the range of effective radiative forcings for all aerosol-radiation interactions is -0.95 to 0.05 W/m 2 5 and for aerosol-cloud interactions is -1.2 to 0.0 W/m 2 . The level of scientific understanding in that report for the cloud albedo effect is still marked as "low". From the air quality perspective, in extreme cases sulfur chemistry may lead to creation of acid rain or acid fog (Dianwu et al., 1988;Wang et al., 2016). Based on analysis of 20 modeling studies, the review by Faloona (2009) marks wet deposition of aerosol sulfate, dry deposition of SO 2 and heterogeneous (aqueous phase) oxidation of SO 2 in aerosol particles and clouds as the most challenging to quantify in models. For an overview of the representation of sulfur 10 oxidation in regional and global models see Ervens (2015). The aqueous phase oxidation is reported as a dominant mechanism of production of sulfate: A numerical study by Barth et al. (2000) reports that for the in-cloud conditions, aqueous phase reactions accounts for 81% of sulfate production rate. According to their study total of ∼ 50% − 60% of sulfate burden in the troposphere is produced by aqueous phase chemistry. The gas phase SO 2 is oxidized in a matter of days by the gas phase reactions or within minutes or a few hours within clouds by aqueous phase reactions, see the review by Faloona (2009). 15 From the cloud microphysics stand point, aqueous phase oxidation of sulfur is interesting because it affects the CCN within water drops. Noteworthy, sulfate is a common component of aerosol particles (10%-67% of the submicron particle mass is made of sulfate, 32% on average; see Zhang et al., 2007). The aqueous phase oxidation of sulfur is irreversible, meaning that the produced sulfate remains within the water drops and increases the dissolved CCN mass. Collisions and the subsequent coalescence of water drops as well as collisions between the aerosol particles and water drops are another in-cloud irreversible 20 process that affects the aerosol particles. As the water drops collide and coalesce, the created new water drop carries within the combined CCN mass of all of its colliding predecessors. Efficient collisions between cloud droplets may quickly lead to the onset of precipitation which can in turn effectively cleanse the atmosphere from aerosol particles and water soluble trace gases.
In non-precipitating clouds, aerosol particles that served as CCN are altered by cloud microphysical and chemical processes and then returned to the atmosphere after water drops evaporate (the process is referred to as the CCN deactivation, aerosol 25 regeneration, aerosol recycling or aerosol resuspension). The cloud-processed aerosol particles can be observed in measurements (Hoppel et al. 1986, 1994, Werner et al. 2014, Hudson et al. 2015. The case without precipitation reaching the surface is especially interesting as it allows for aerosol-cloud interactions to loop for several cloud life-cycles without removing the altered aerosol particles. The 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 30 aerosol particle has been cycled 3 times by cloud systems. The Particle-Based Microphysics and Chemistry (PBMC) scheme introduced here offers a chance to represent the effects of such cloud-processing on CCN sizes stemming from both collisions between water drops and aqueous phase oxidation reaction within water drops. The PBMC can be used in multidimensional simulations with a fully coupled dynamics model, which has not been possible before. To the authors knowledge, the presented 2 Intergovernmental Panel on Climate Change, see http://www.ipcc.ch/ scheme is the first to represent the impact of both collisions and aqueous phase chemistry on the aerosol size spectrum in the particle-based microphysics framework. This paper documents the extension of the particle-based microphysics scheme with a numerical scheme that represents aqueous phase chemical reactions inside cloud droplets and the uptake of the trace gases into cloud droplets. The representation 5 of chemical reactions includes only the aqueous phase processes (i.e., no gas phase chemical reactions) and revolves around oxidation of sulfur dissolved in water drops to sulfate. Two reaction paths are considered -oxidation by ozone and by hydrogen peroxide. In total, six trace gases are included in the chemistry description: sulfur dioxide (SO 2 ), ozone (O 3 ), hydrogen peroxide (H 2 O 2 ), carbon dioxide (CO 2 ), nitric acid (HNO 3 ), and ammonia (NH 3 ). Their dissolution, and if applicable, dissociation is resolved. The structure of the presented work is as follows: Section 2 presents briefly the particle-based scheme available in 10 libcloudph++. Section 3 discusses the design and programming 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 PBMC 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.
2 Particle-based microphysics scheme 15 The particle-based 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 particle-based 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 super-droplets (SD). Each SD has a set of attributes describing the properties of the aerosol particles or water drops it 20 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 3 (r d ) and the hygroscopicity parameter 4 (κ).
The aqueous phase chemistry scheme extends the list of required attributes by masses of chemical compounds dissolved in droplets. Eight new attributes are needed, see Sec. 3 for details.
The key attribute of the particle-based microphysics scheme is the SD multiplicity. Multiplicity defines the number of aerosol 25 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 particle-based 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 30 of condensational growth is solved for each SDs wet radius (see Sec. 5.1.3 in Arabas et al. (2015) for details). The process of 3 It is a volume equivalent radius for solute in the water drop. 4 Following Ghan et al. (2001) and Petters and Kreidenweis (2007) it is a single parameter representing the hygroscopicity of the solvent. In this work we use the notation from Petters and Kreidenweis (2007). 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, see 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.4 in 5 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 particle-based 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 particle-based approach 10 greatly reduces the numerical diffusion errors. As discussed in Unterstrasser et al. (2017), 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. Dziekan and Pawlowska (2017) showed that for high SD concentrations the super-droplet method accurately represents collisions between the drops (with regard to the expected value and the standard deviation of 15 the autoconversion time). An interesting comparison between the bin and the particle-based schemes is provided by Li et al. (2017). On a side note, Grabowski and Abade (2017) show that an additional scheme modeling the broadening of droplet spectra due to supersaturation fluctuations might be necessary.
The collision efficiency used in this study is based on Hall (1980) and Pinsky et al. (2008). It is well suited for representing the collisions between water drops. An additional collision efficiency look-up table based for example on Ladino et al. (2011) 20 or Ardon-Dryer et al. (2015) should be used to study the collection of submicron aerosol particles by droplets. Similarly, additional collision efficiency corrections based for example on Chen et al. (2018) should be applied to study the effects of turbulence on the aerosol size distribution.
3 Aqueous phase chemistry scheme 30 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 acidity. This results in eight new SD attributes needed for simulations with aqueous phase chemistry: the total mass of dissolved O 3 , the total mass of dissolved H 2 O 2 , the total mass of dissolved SO 2 (including: SO 2 * H 2 O, HSO − 3 and SO 2− 3 ),

10
The scheme needs to be coupled to a driver model (i.e. dynamical core) that provides information about the environment in which SD are immersed (i.e. temperature, humidity, trace gas mixing ratio, and air wind field). 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 15 any additional parametrisation.
All aqueous-phase chemistry included in the scheme is formulated under the assumption that solution droplets are diluted. Therefore, in the PBMC 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 20 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 particle-based 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 recalcu- 25 lating 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.

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): 5 where ∆H D denotes the reaction enthalpy of dissociation at constant temperature and pressure, T is the temperature of air and R is the gas constant. The list of considered dissociation constants as well as their temperature dependence coefficients is available in Tab. C2. The dissociation of water, although very small, is also taken into account 5 . It is assumed that the water dissociation constant does not vary with temperature. 10 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 . (2) 15 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 the respective 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 20 obtained iteratively using a numerical root-finding algorithm 6 . 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. 25 When the SD wet radius is quickly changing, for example during the initial condensational growth of cloud droplet or rain 5 The concentration of undissociated water molecules is so big that it is usually assumed constant and it traditionally multiplies the dissociation constant of water. This leads to a different definition of the dissociation constant for water: H2O = [H + ][OH − ] 6 TOMS 748 algorithm from Boost library. See www.boost.org/doc for documentation and Alefeld et al. (1995) for derivation 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 5 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 where D A and α M A are the diffusion and mass accommodation coefficients for the chemical compound "A", <v> = 8RT  Table C3 shows the Henry's law constants and their temperature dependencies 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 (2016) chap. 7.3 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 20 ambient trace gas concentration is calculated from the trace gas mixing ratio provided by the driver model to which the PBMC 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 PBMC 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. To ensure that the sum of sinks from each SD does not exceed the available ambient trace gas 25 mixing ratio, a 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 dependence coefficients.

5
Equations (4) and (5) return the new concentration of S VI created in each SD in each time-step. Based on the new concentration, the new mass of S VI 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 S VI mass.
For the typical atmospheric conditions, say pH between 3 and 6 (i.e. [H + ] between 10 −3 M and 10 −6 M), it can be said that

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 15 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 set to zero and is therefore not in equilibrium with the initial ambient trace gas conditions. For the initial conditions where supersaturation is present in the environment 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. Such initial conditions are mostly relevant for the kinematic models. 20

Validation
The PBMC scheme is set to reproduce results from model intercomparison study by Kreidenweis et al. (2003), where several bulk and moving-bin schemes representing cloud microphysics and aqueous-phase chemistry were tested in an adiabatic parcel model setup. The parcel model used here 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 super- 25 saturation. 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 30 starts below cloud base (i.e. with subsaturation). The initial aerosol is ammonium bisulfate and the initial aerosol particle size 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. See Sec. 5.1.6 in Arabas et al. (2015) for the details on how the super-droplet dry and wet radii are initialized. 5 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 needs to assume a vertical profile of pressure. In the presented simulations the pressure profile is obtained by integrating the hydrostatic equation and assuming that density of air is constant and equal to 1.15 kg m -3 . The assumed density is based on the density provided in Tab. 3 in Kreidenweis et al. (2003). Then, at each level, ρ d is calculated from the ideal gas law taking into account the current r v and θ: where: p v represents partial pressure of water vapor, p 1000 stands for pressure equal 1000 hPa that comes from the definition of potential temperature, R d is the gas constant for dry air and c pd is the specific heat at constant pressure for dry air. Because the simulated air parcel is assumed to be adiabatic, only processes resolved by the particle-based scheme can change θ, r v and 5 other trace gas mixing ratios. In each model time-step, the particle-based 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 10 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 volume 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 15 compares well with Fig. 1 from Kreidenweis et al. (2003).
At the end of the test simulation, 85% of SO 2 is converted into S VI and the final water volume weighted average pH is equal to 4.86. The total sulfate production is 171 ppt with 99 ppt produced by the H 2 O 2 reaction path and 72 ppt produced by the O 3 reaction path. Based on Fig. 2 in Kreidenweis et al. (2003), the range of average pH values reported by different size resolving (moving-bin) schemes was between 4.82 and 4.85, and the range of total sulfate production values was between 170 -180 ppt.
Based on Fig. 3 in Kreidenweis et al. (2003) the production by H 2 O 2 ranged between 85 and 105 ppt, and by O 3 between 70 and 85 ppt for the size resolving schemes. In short, the results from the particle-based scheme are close to the range of values reported by moving-bin schemes. 5 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. 2 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 moving-bin schemes responsible for causing this discrepancy is the different water vapor mass accommodation coefficient α M leading to different predicted 10 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 particle-based scheme used in this study reports concentration 269 cm −3 at the level of maximum supersaturation. The maximum supersaturation is equal to 0.27%. The particle-based scheme assumes α M equal to unity and therefore fits with the trend of high α M causing lower supersaturation presented in Fig. 6 in 15 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, see also the discussion in Arabas and Pawlowska (2011). Along those lines, here it is tested how sensitive the particle-based 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  particles grow in size due to S VI production during oxidation, but the increase in size is small compared to their initial size. The smallest activated aerosol particles are those that are affected most by oxidation. The increase in their size due to the produced S VI is the largest compared to their initial size. In short, oxidation produces a "gap", often labeled the "Hoppel minimum", between the CCN processed by the cloud and the smaller unactivated aerosol particles.
The effect of in-cloud sulfate production on aerosol particle size distribution presented in Fig. 3, combined with other tests 5 presented in this chapter, document the correctness of the implementation of aqueous chemistry in the particle-based scheme.
The formation of the "Hoppel minimum" was reported by many observational studies, see Hoppel et al. (1994), Bower et al. (1997) and 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), Kreidenweis (2000, 2002) and Ovchinnikov and Easter (2010). The work by Cantrell et al. (1999) shows that the maximum supersaturation and the cloud droplet concentration in the clouds which processed the aerosol particles can be inferred based on the location of the aerosol size distribution minimum. The collisions between water drops are represented using the geometric kernel with collision efficiency for big drops (i.e. 15 radius greater than 20 µm) from Hall (1980) and for small droplets (i.e. radius smaller than 20 µm) 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 5 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 10 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 model variables. 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 20 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), sea salt aerosol particles are alkaline, which may in turn increase the pH of water drops and thus affect the oxidation rate. On 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 35  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. unactivated aerosol concentration (a), cloud droplet specific concentration (b), rain water mixing ratio (c), mean dry radius (d), cloud droplet effective radius (e) and concentration of S VI molecules (f). The thresholds for particle radii are: unactivated aerosol < 1 µm; 1 µm < cloud < 25 µm; rain > 25 µm. Note the logarithmic scale for the rain water mixing ratio plot.
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 5 straightforward selection of the updraft and downdraft regions, further simplifying the analysis of the microphysical processes. 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 in-cloud aerosol particles that did not activate. The difference between the 5 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 10 regions with drizzle (see Fig. 4c). 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. In the updraft (left-hand 15 side of panel d) environmental aerosol particles that have not been affected by cloud are advected into the cloudy region. Once the cloud droplets are formed, the aqueous phase oxidation starts to produce sulfate and changes the CCN size distribution. In the downdraft (right-hand side of panel d) cloud droplets are advected out of the cloud and they evaporate. The cloud-processed CCN are returned to the environment and change the ambient air aerosol 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 20 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)). 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 25 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 grid-cells. Figure 4f also shows that the particle-based scheme can track the dissolved chemical compounds in the evaporating rain drops below the cloud base. 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 30 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.

Results
The depicted clean layer is an artifact caused by the prescribed velocity field and the absence of aerosol sources in the compu-tational domain. The relatively short simulation time is chosen to minimize the impact of the clean layer on the simulation.  Figure 5 shows the liquid water volume 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 acidic aerosol particles below the cloud base have been clipped. Figure 5 captures the pH of cloud droplets as well 5 as the pH of some evaporating rain drops below the cloud base. The droplets in the downdraft of Fig.5a have higher acidity that is caused by S VI created during aqueous phase oxidation. For base case, Fig. 5a, pH increases with height above the cloud base. Initially water drops are very acidic, but they grow in size and become more diluted. Even though S VI is created during oxidation, the average pH still increases with height due to dilution. The same behavior is shown in the adiabatic parcel tests discussed in Sec. 4 and shown in Fig. 1. The increase of pH with height is also observed in 1-dimensional model representing 10 processing of sulfur in small cumuli in marine environment (Alfonso and Raga, 2002). Due to the pH variability shown in increases significantly with increasing pH (see Eq. 4) whereas oxidation by H 2 O 2 depends very weakly on the acidity (Eq. 5).
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 5 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 10 volume 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, see Eq.
(2). Then, as the water drops are advected upwards, oxidation produces sulfuric acid and the average pH decreases. Near the cloud top, the 15 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 20 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 25 minutes simulation time (not shown). Figure 5d shows the liquid water volume 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 less diluted. 30 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 S VI originates from S IV oxidation by H2O2. 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 5 case3 simulation, 39% of the S VI originates from S VI oxidation by O 3 and 61% by H 2 O 2 . 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 10 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). 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 15 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 bi-modal 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 5 because the produced sulfate is divided among bigger number of aerosol particles.

Summary and outlook
The work presented here describes an extension of the libcloudph++ library that allows including aqueous phase chemical reactions within water drops in the particle-based 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, 10 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 validate the the particle-based microphysics and chemistry scheme. Additionally, the changes in the programming 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. 15 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 particle-based scheme but do not provide a good balance between the representation of cloud microphysics and dynamics. As a next step, the particle-based 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. 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  plemented 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 particle-based microphysics scheme in lib-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 particle-based microphysics scheme options are grouped into a struc-5 ture 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 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 10 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 in Fig. A2 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) arrinfo_t<real_t> rv, const arrinfo_t<real_t> rho const arrinfo_t<real_t> cou const arrinfo_t<real_t> cou const arrinfo_t<real_t> cou std::map< enum chem_species_t, arrinfo_t<real_t> > ambient_chem, ); void step_async( const opts_t<real_t> & ); // diagnostic methods // ... void diag_chem(const enum che // ... Listing 3.4: lgrngn::particl 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 8 Figure A3. lgrngn::chem_species_t definition The names of chemical compounds available in the aqueous phase chemistry module are stored in a chem_species_t enum, see code listing in Fig. A3. The state of all variables used by the particle-based 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 CPU and GPU calculations (see the discussion in Sec. 5.2 in Arabas et al., 2015, for details). The initialization, time-stepping and output from the particle-based scheme are done using the methods of lgrngn::particles_t structure. Their 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. 10 They are optional and depend on the dimensionality of the solved problem. For example, for the parcel model tests from Sec. 4 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 enum 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. 15 During time-stepping the particle-based scheme computations are performed by step_sync() and step_async() methods.
The first one gathers all the processes that affect the driver model fields (such as condensation/evaporation or aqueous phase chemistry) 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 20 map linking chem_species_t enum items with the driver model mixing ratio fields needs to be provided as the last optional argument. The particle-based scheme overwrites the driver model fields during simulation. The signature of step_async method is not changed by the new chemistry module. Listing 3.4: lgrngn::particles_t definition 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 As discussed in Arabas et al. (2015), the lgrngn::particles_t structure provides many methods for obtaining statistical information on the SD parameters (prefixed with diag). The chemistry model adds to them the diag_chem method that outputs the total mass of a chemical compound dissolved into droplets. The chemical compound is selected using the chem_species_t enum items. See the discussion in Sec. 5.2 in Arabas et al. (2015) for the details on how to select the size ranges of droplets specified for output or how to output other statistical parameters. 5  Table C2. Dissociation constants and their temperature dependence coefficients (taken from Kreidenweis et al., 2003). Dissociation of H 2 SO 4 is taken from Tab. 7A.1 in Seinfeld and Pandis (2016). equilibrium reaction dissociation constant at 298K (moles liter -1 ) temp. dep. −∆H D