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.

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. particle tracking or “Lagrangian-in-droplet-radius-and-space”) microphysics scheme.
The scheme resolves the evolution
of the aerosol, cloud droplet, and rain drop1
size spectrum.
It allows for representing cloud microphysical processes from first principles
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-source and its programming interface
is documented in Arabas et al. (2015).
All those features facilitate further development and usage
of libcloudph++.

Sulfate aerosols cool Earth's climate by scattering sunlight – and thus increasing
Earth's shortwave albedo (direct radiative forcing) – and also by changing radiative properties of clouds
(cloud albedo effect).
According to chapter 8 of IPCC2
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 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, sulfur chemistry may lead to the creation of
acid rain or acid fog in extreme cases (Dianwu et al., 1988; Wang et al., 2016).
Based on analyses of 20 modelling studies,
the review by Faloona (2009) marks wet deposition of aerosol sulfate,
dry deposition of SO2, and heterogeneous (aqueous phase)
oxidation of SO2 in aerosol particles and clouds
as the most challenging to quantify in models.
For an overview of the representation of sulfur oxidation in
regional and global models see Ervens (2015).
Aqueous-phase oxidation is reported as a dominant mechanism of the production of sulfate:
a numerical study by Barth et al. (2000) reports that for the
in-cloud conditions, aqueous-phase reactions account for 81 % of the sulfate production rate.
According to their study, a total of ∼50 %–60 % of sulfate burden in the troposphere
is produced by aqueous-phase chemistry.
The gas-phase SO2 is oxidized in a matter of days by gas-phase
reactions, or within minutes or a few hours within clouds by aqueous-phase reactions,
see the review by Faloona (2009).

From the cloud microphysics stand point, aqueous-phase oxidation of sulfur is interesting
because it affects the CCN within water drops.
Note that sulfate is a common component of aerosol particles
(Zhang et al., 2007).
The aqueous-phase oxidation of sulfur is irreversible, meaning that the produced sulfate
remains within water drops and increases the dissolved CCN mass.
Collisions and the subsequent coalescence of water drops
as well as collisions between aerosol particles and water drops
are another in-cloud irreversible process that affects aerosol particles.
As water drops collide and coalesce, the newly created water drop
carries 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
return to the atmosphere after water drops evaporate
(the process is referred to as the CCN deactivation, aerosol regeneration, aerosol recycling, or aerosol resuspension).
The cloud-processed aerosol particles can be observed in measurements
(Hoppel et al., 1986, 1994; Hudson et al., 2015; Werner et al., 2014).
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 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 reactions 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 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 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 – the oxidation by ozone and by hydrogen peroxide.
In total, six trace gases are included in the chemistry description:
sulfur dioxide (SO2), ozone (O3), hydrogen peroxide (H2O2),
carbon dioxide (CO2), nitric acid (HNO3), and ammonia (NH3).
Their dissolution and, if applicable, dissociation is resolved.
The structure of the presented work is as follows:
Sect. 2 briefly presents
the particle-based scheme available in libcloudph++.
Section 3 discusses the design
of the new aqueous chemistry scheme and Sect. A in the Appendix
describes the new programming interface.
Section 4 compares the results from the new scheme
with the results from moving-bin schemes.
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.

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 modelling cloud microphysics,
the computational domain is filled with “numerical point particles”
representing a specified number (called 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 labelled here as super-droplets (SDs).
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
the multiplicity (𝒩),
the position of SD in the computational domain,
the wet radius (rw),
the dry radius3
(rd),
and the hygroscopicity parameter4
(κ).
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 Sect. 3 for details.

The key attribute of the particle-based microphysics scheme is the SD multiplicity.
The 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.
The use of multiplicity reduces 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, or rain water,
as it is often done in bulk schemes,
for example Kessler (1995), Seifert and Beheng (2001), and Morrison and Grabowski (2007).
All the modelled microphysical processes are represented by calculating
the changes to the SD attributes.
The equation of condensational growth is solved for each SDs wet radius
(see Sect. 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 parameterization of cloud droplet activation is required
as it is again often done in bulk microphysics schemes, see for example
Morrison and Grabowski (2007).

Following Shima et al. (2009),
the collisions between SDs are represented using a Monte Carlo scheme
(see Sect. 5.1.4 in Arabas et al., 2015, for details).
Information about the 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 rw, rd, and 𝒩.
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 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 SD method accurately represents collisions between the drops
(with regard to the expected value and the standard deviation of 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 modelling the broadening of droplet spectra
due to supersaturation fluctuations might be necessary
for the particle-based schemes.

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 on, for example,
Ladino et al. (2011) 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 on, for example,
Chen et al. (2018) should be applied
to study the effects of turbulence on the aerosol size distribution.

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 the SD's acidity.
This results in eight new SD attributes needed for simulations with aqueous-phase chemistry:

the total mass of dissolved O3,

the total mass of dissolved H2O2,

the total mass of dissolved SO2 (including: SO2*H2O, HSO3-, and SO32-),

the total mass of dissolved CO2 (including: CO2*H2O, HCO3-, and CO32-),

the total mass of dissolved NH3 (including: NH3*H2O and NH4+),

the total mass of dissolved HNO3 (including: HNO3(aq) and NO3-),

the total mass of created H2SO4 (including: HSO4- and
SO42-), and

the total mass of H+ ions.

The scheme needs to be coupled to a driver model (i.e. a dynamical core) that
provides information about the environment in which SDs 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 any additional parameterization.

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 the SDs
with ionic strength smaller than 0.02 moles L−1
(Ovchinnikov and Easter, 2010).
In practice, this condition results in excluding
the SDs with small wet radii from aqueous chemistry calculations
(i.e. the SDs representing haze particles and very small cloud droplets).
Exclusion of the SDs with small wet radii also
prevents numerical issues during the 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 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.

3.1 Dissociation

Dissociation is a reversible process of splitting 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):

(1)KA(T)=KA(T0)exp-ΔHdR1T-1T0,

where ΔHD 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
and their temperature dependence coefficients are available in Table C2.
The dissociation of water, although very small, is also taken into account5.
It is assumed that the water dissociation constant
does not vary with temperature.

It is assumed that there is no electric charge in water drops
and therefore the concentrations of positive and negative ions
created during dissociation should balance each other.
Using the dissociation constants (see Table 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.
For an example derivation, see Sect. 7.6.2 in Seinfeld and Pandis (2016).
The neutral charge condition can be expressed as

The [] brackets denote the concentration of each of the chemical compounds
(traditionally defined in units of moles per litre),
capital letters denote the chemical compound and
roman numbers mark its oxidation state.
In Eq. (2) the dissociation constants of SO32-, CO32-, and SO42- ions
(i.e. KHSO3, KHCO3, and KH2SO4)
are multiplied by a factor of 2 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 obtained iteratively
using a numerical root-finding algorithm6.
The algorithm searches for a solution between pH=-1 and pH=9.
The lower bound for the pH scale is unrealistically low
and is only necessary at the start of the simulation
when the initial SDs have very small volume and are highly acidic.
The upper bound is set arbitrarily, but is sufficient
for the expected pH of the modelled 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 the dissociation procedure can be divided into a user-specified number of sub-steps
in order to prevent limiting the overall simulation time step by dissociation.

3.2 Dissolution

The amount of chemical compound that can dissolve into a 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 phases
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 is given
by

(3)d[A]dt=4rw3〈v〉αMA+rw23DA-1cA∞-[A]HAeffRT

where DA and αMA are
the diffusion and mass accommodation coefficients for the chemical compound
“A”, 〈v〉 = 8RTπMA is the average velocity of the molecules
calculated from the Maxwell–Boltzmann distribution function,
MA is the molar mass
of the chemical compound “A”, cA∞ is the ambient concentration of the trace gas “A”, and
HAeff is the effective Henry's
law constant of the chemical compound “A” (i.e. the equilibrium dissolution constant).
The Henry's law constants depend on the temperature following a similar relation as
for dissociation, Eq. (1).
Table C3 shows the Henry's law constants
and their temperature dependencies and
Table C4 presents the diffusion and mass accommodation coefficients.
The term “effective” marks that the dissolution constants take into account
the increase in the efficiency due to dissociation
(see Seinfeld and Pandis, 2016, Sect. 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 the 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 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 the 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 mixing ratio,
a relatively short time step should be applied.
If necessary, the user can divide the model time step into sub-steps.

where ℝA is the reaction rate of the chemical
compound “A” and k0,…,4 are the reaction rate coefficients.
k0,…,4 depend on the temperature following a similar relation as
for dissociation, Eq. (1).
Table C5 shows the values of the reaction rate coefficients and their
temperature dependence coefficients.

Equations (4) and (5)
return the new concentration of SVI created in each SD in each time step.
Based on the new concentration,
the new mass of SVI
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 SVI mass.

For the typical atmospheric conditions, say pH between
3 and 6 (i.e. [H+] between 10−3M and 10−6M),
it can be said that the rate of oxidation by H2O2 depends very weakly on pH.
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 unit results in an approximately 100-fold increase in the O3
reaction rate.

3.4 Initialization

The initial aerosol is assumed to be ammonium bisulfate (NH4HSO4),
with a dry particle density of 1.8 g cm−3.
Using the dry particle density and the dry radius of each SD,
the initial mass of H+, NH4+ and SO42- 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.

The PBMC scheme is set to
reproduce results from the 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 set-up.
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 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 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

(6)n(rd)=ntotrd2πln(σg)exp-ln(rd)-ln(rd‾)22ln2(σg),

where n(rd) is the spectral density function of aerosol particle sizes,
ntot is the total aerosol concentration,
rd‾ is median radius, and
σg is the geometric standard deviation.
See Sect. 5.1.6 in Arabas et al. (2015) for the details on how the
SD dry and wet radii are initialized.

The parcel model employed in this study uses dry air density
ρd, dry air potential temperature θ, water vapour mixing ratio rv,
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 the density of air is constant and equal to
1.15 kg m−3.
The assumed density is based on the density provided in Table 3 in Kreidenweis et al. (2003).
Then, at each level, ρd is calculated from the ideal gas law
taking into account the current rv and θ:

(7)ρd(p,θ,rv)=p-pv(p,rv)Rdθ(pp1000)Rdcpd,

where:
pv represents the partial pressure of water vapour,
p1000 stands for the pressure equal 1000 hPa that comes from the definition of potential temperature,
Rd is the gas constant for dry air, and
cpd is the specific heat at constant pressure for dry
air. Because the simulated air parcel is assumed to be adiabatic,
only the processes resolved by the particle-based scheme can change θ,
rv, and other trace gas mixing ratios.
In each model time step, the particle-based microphysics scheme changes θ and rv
according to Eqs. (25) and (26) from Arabas et al. (2015).
The changes in the trace gas mixing ratios are resolved
following the procedure discussed in Sect. 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−1 at a height of 1.2 km above the cloud base.
Figure 1b shows the total SO2 concentration
(both in gas phase and dissolved in water) in ppb units.
The concentration of SO2 is decreasing due to oxidation
taking place in the cloud droplets.
Figure 1c shows the water volume weighted average pH
of the cloud droplets.
The pH near the cloud base is very low due to the acidic nature of the assumed
initial aerosol and the small size of the 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).

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

At the end of the test simulation,
85 % of SO2 is converted into SVI
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 H2O2 reaction path
and 72 ppt produced by the O3 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 and 180 ppt.
Based on Fig. 3 in Kreidenweis et al. (2003),
the production by H2O2 ranged between 85 and 105 ppt,
and by O3 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 the moving-bin schemes.

The microphysics schemes taking part in the Kreidenweis et al. (2003) intercomparison study reported
significant differences between the number of activated cloud droplets.
Based on Table 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 vapour 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 % and 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 a concentration of
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 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 total sulfate production (d).
The results are plotted against the logarithm of base 2 of the number of SDs
in the computational domain (meaning that “0” represents one SD
and “10” represents 1024 SDs).
All values seem to converge for SD numbers greater than 64.
The average pH, maximum supersaturation, and total sulfate production
do not change for those four test-runs.
The concentration of droplets at the cloud base varies little
(between 269 and 281 cm−3).
The concentrations from simulations with SD number between 512 and 1024
vary between 274 and 269 cm−3.
For SD numbers between 32 and 64 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 281 and 332 cm−3).
This confirms the observations from Kreidenweis et al. (2003)
that the predicted cloud droplet number concentration
strongly depends on the representation of the size distribution of
modelled 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 2
Results of the convergence test for the adiabatic parcel simulations.
All panels show how a given parameter depends on the number
of SDs (shown on the abscissa as the logarithm of base 2 of the number of SDs).
Panel (a) shows the cloud droplet concentration at the cloud base,
panel (b) the maximum supersaturation,
panel (c) the water volume weighted average pH at the end of
simulation, and panel (d) the total sulfate production.

Figure 3 shows
the simulated modification of the aerosol size distribution.
The red line depicts the initial distribution
and the green line shows model state at the end of the adiabatic parcel test.
For convenience, Fig. 3 uses both logarithmic (left panel) and linear (right panel)
scales on the axes.
The change in the aerosol size distribution is caused by oxidation.
The aerosol particles that are too small to become cloud droplets
are not affected by aqueous-phase oxidation and do not grow in size.
The large aerosol particles grow in size due to SVI production during oxidation,
but the increase in size is small compared to their initial size.
The smallest activated aerosol particles
are affected most by oxidation.
The increase in their size due to the produced SVI is the largest
compared to their initial size.
In short, oxidation produces a “gap”,
often labelled the “Hoppel minimum”,
between the CCN processed by the cloud and the smaller
unactivated aerosol particles.

Figure 3
Modification of the dry aerosol sulfate mass.
The red line shows the initial condition and the green line shows the final model state.
The left panel uses a logarithmic scale and the right panel uses a linear scale.

The effect of in-cloud sulfate production on the aerosol particle size distribution
presented in Fig. 3,
combined with other tests presented in this section,
documents the correctness of the implementation of the 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 the 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 the Hoppel minimum,
see for example Flossmann (1994),
Feingold and 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.

5.1 2-D kinematic model

The kinematic model mimics a single 2-D eddy
spanning a stratocumulus cloud deck and a boundary layer below it.
The model is based on a test scenario from
the 8th International Cloud Modeling Workshop (Muhlbauer et al., 2013).
The velocity field is prescribed as in Szumowski et al. (1998),
Morrison and Grabowski (2007), and Rasinski et al. (2011).
The same model was used when presenting the initial release of libcloudph++,
see Sect. 2 in Arabas et al. (2015) for details of the model formulation.
The kinematic model is based on the open-source library of parallel MPDATA-based
solvers for systems of generalized transport equations, see Jaruga et al. (2015).
The temperature, moisture, and trace gas fields are
discretized on the Eulerian grid and are
advected using the prescribed velocity field.
Then, the model variables are passed to the PBMC scheme,
where the microphysical and chemical processes are resolved.
Finally, the source and sink terms due to microphysics and chemistry
are calculated and applied in each model grid cell as described in Sects. 2 and
3.

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

The collisions between water drops are represented using the
geometric kernel with collision efficiency
for big drops (i.e. 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 table for different drop sizes.

The initial conditions are summarized in Table 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 rv are constant with altitude.
To keep the simulation set-up 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 105 s (∼2 h 45 min) 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 % (relative
humidity RH=1.05).
After spin-up the simulations are run for 30 min.
The chosen simulation time is enough to deplete the SO2 available in the cloudy
part of the domain as well as to create precipitation.

Similar 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 initial SO2, O3, and H2O2 volume fractions are
taken from the simulation set-up used in Ovchinnikov and Easter (2010).
The values for SO2 and O3 are based on the
measurements from the MASE campaign (Wang et al., 2008)
and the value for H2O2 is based on the representative values
for the eastern Pacific Ocean (Genfa et al., 1999).
The NH3, HNO3, and CO2 volume fractions
are the same as in the parcel test from Sect. 4.

The set-up detailed in Table 2 corresponds to “very clean conditions”
(i.e. low aerosol particle concentrations).
The 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 the oxidation reaction is prohibited.
In case2 the initial volume fraction of NH3 is increased and in
case3 the initial aerosol size distribution is changed.
The conditions for all the sensitivity simulation cases
are summarized in Table 3.

Table 3Initial conditions for sensitivity test cases of the 2-dimensional kinematic model.
Specified are aqueous-phase chemistry choice,
initial volume fraction of NH3, mean radius of the assumed
lognormal aerosol particle size distribution rd, total aerosol concentration
ntot,
and geometric standard deviation σg.
Other parameters for each case are the same as in the base case (Table 2)
The parameters that distinguish each sensitivity test case
are marked in bold.

As discussed in Flossmann (1994),
the initial chemical scenario is idealized.
For instance, although the initial conditions represent a clean maritime environment,
the set-up 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 other hand, a study by von Glasow and Sander (2001)
indicates that alkaline sea salt particles are quickly converted to acidic
due to the uptake of HCl vapour.
More importantly, including sea salt would result in
aerosol particles with very different hygroscopicity values
(Petters and Kreidenweis, 2007).
Including sea salt would also result in the initial bimodal size distribution
with one mode representing smaller ammonium bisulfate aerosol particles
and the second mode representing larger sea salt particles.
In general, including sea salt should result in a very different condensational
growth of aerosol particles.
The set-up used in this study also lacks other particles containing sulfate,
such as ammonium sulfate or sulfuric acid aerosol particles.

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 rd‾ and precipitation:
(i) for a relatively small initial rd‾
production of sulfate enhances precipitation, and
(ii) for a relatively big initial rd‾
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 large eddy simulations
should focus on the impacts of aqueous chemistry on precipitation, cloud lifetime, and cloud dynamics.

The kinematic set-up precludes any links between cloud microphysical processes and dynamics of the air motion.
The set-up 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 set-up has low computational cost
and allows for 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 set-up allows for a straightforward selection of the updraught
and downdraught regions, further simplifying the
analysis of the
microphysical processes.

Figure 4 shows the model state after 30 min of simulation from the base case
(see Table 2).
Figure 4a shows the concentration of the 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 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 the cloud droplets.
Figure 4b shows the concentration of the 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
(Pawlowska et al., 2000).
The regions with lower cloud droplet concentrations correspond
to the 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 colour 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. 4b and f
where the cloud droplet concentration is depleted in regions of drizzle.
Figure 4d shows
the mean dry radius of all particles (both the aerosol particles and water drops).
The mean dry radius is increasing due to oxidation.
In the updraught (left-hand side of panel d) the 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 downdraught (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 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 Pawlowska and Brenguier, 2003; Rosenfeld and Gutman, 1994).
Figure 4f shows the concentration of SVI molecules
(all molecules containing sulfur at
+VI oxidation state)
and represents molecules from the initial ammonium bisulfate (NH4HSO4) 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 SVI 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 the downdraught 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.
The depicted clean layer is an artifact caused by the prescribed velocity
field and the absence of aerosol sources in the computational domain.
The relatively short simulation time
is chosen to minimize the impact of the clean layer
on the simulation.

Figure 5
The liquid water volume weighted average pH from
the base case (a),
case1 (b), case2 (c), and case3 (d).
See Tables 2 and 3 for a definition of simulation set-ups.

Figure 5 shows the liquid water volume weighted average pH
in each computational grid cell from the base case (a) and sensitivity test cases (b–d).
In order to better adjust the colour 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 as
the pH of some evaporating rain drops below the cloud base.
The droplets in the downdraught of Fig. 5a
are more acidic due to SVI created
during aqueous-phase oxidation.
For the base case, Fig. 5a, pH increases with height above the cloud base.
Initially the water drops are very acidic, but as they grow in size they
become more diluted.
Even though SVI is created during oxidation,
the average pH still increases with height due to dilution.
The same behaviour is shown in the adiabatic parcel tests
discussed in Sect. 4 and shown in Fig. 1.
The increase in pH with height is also observed in a
1-dimensional model representing processing of sulfur
in small cumuli in marine environment used by Alfonso and Raga (2002).
Due to the pH variability shown in Fig. 5a,
oxidation by O3 happens mostly near the cloud top
in the base case.
As discussed in Sect. 3.3, the rate of oxidation by O3
increases significantly with increasing pH (see Eq. 4),
whereas oxidation by H2O2 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 O3 in the upper parts of the cloud,
whereas oxidation by H2O2 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 the base case,
the reversible chemical processes are allowed, and the oxidation reaction is prohibited.
The scenario without oxidation
is overall less acidic than the base case (Fig. 5a).
Additionally, without oxidation there is no difference between the pH values
in the updraught and downdraught 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 downdraught.

Figure 6
The size distributions of dry radii for the base case (a) and case3 (b).
The initial dry radius size distribution is marked in black, the final dry radius size distribution
from grid cells with rc>0.01 g kg−1 in green, and from
grid cells with rr>0.01 g kg−1 in red.
See Tables 2 and 3 for a
definition of simulation set-ups.

Case2 differs from the base case by increasing the initial
NH3 volume fraction from 0.1 to
0.4 ppbv
(see Tables 2 and 3).
Because the initial aerosol particle size distribution is
the same as in the base case,
the mean aerosol and droplet sizes and concentrations at the end
of the simulation are not different from the base case (not shown).
Figure 5c shows the liquid water volume weighted average pH for case2.
The average pH in case2 (Fig. 5c) is higher than in the base case (Fig. 5a),
that is, both cloud droplets and rain drops are less acidic in case2 than in the base case.
In contrast to the base case, in the updraught (left-hand side of the plots),
the pH in case2 actually decreases with height above the cloud base.
This is because the higher initial NH3 volume fraction increases its uptake and
counters the low pH values caused by the 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 NH3 is degassed back to the environment.
The 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 the 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 NH3 volume fraction
is 10 times larger than the base case value assumed here).
As in case2 presented here,
the high initial NH3 volume fractions in Zhang et al. (1999)
increase the pH near the cloud base and promote oxidation by O3
during the first minutes after the simulated parcels entered the cloud.
Because the sulfuric acid was produced, the pH dropped and oxidation by H2O2
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 the base case
(see Tables 2 and 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 min simulation time (not shown).
Figure 5d shows the liquid water volume weighted average pH for case3.
Similar to the base case (Fig. 5a),
the pH increases with height due to the dilution
and the downdraught droplets are more acidic due to the ongoing oxidation.
However, case3 is more acidic than the base case
because the overall droplet sizes are smaller and they are therefore less diluted.

At the end of the base case simulation, 18 % of the total available SIV
is oxidized.
As a result, 0.14 µg m3 of the dry particulate matter
is created during oxidation per cubic metre of air (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 the dry particulate matter
because the initial aerosol mass is very low (only 0.21 µg m3 dry air).
Oxidation by H2O2 is the dominant path:
92 % of SVI originates from SIV oxidation by H2O2.
More alkaline conditions of case2 enhance the efficiency of oxidation.
At the end of case2, 21 % of available SIV is oxidized.
As a result, oxidation produces 0.16 µg of dry particulate matter per cubic metre 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.
Similar to the 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 O3 becomes more
important than in the base case.
At the end of the case3 simulation, 39 % of the SVI originates from SVI oxidation
by O3 and 61 % by H2O2.
In contrast, more acidic conditions of case3 hinder the O3 reaction path.
Virtually all molecules of sulfate that are created during oxidation
are oxidized by H2O2.
As a result, the conversion of sulfur to sulfate is slightly less effective
in case3.
At the end of the case3 simulation, 17 % of available SIV is oxidized.
As a result, 0.13 µg of dry particulate matter is created per cubic metre 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 the base case
due to the higher initial aerosol concentration
(case3 contains initially 0.61 µg m3 of dry particulate matter).
Due to the simple kinematic set-up chosen in this study
the values reported here cannot be treated as representative
of the atmospheric conditions.
They are shown to allow comparison between the 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 the 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 the in-cloud (rc>0.01 g kg−1) and precipitating
(rr>0.01 g kg−1) grid cells, respectively.
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 Sect. 4 and shown in Fig. 3.
Moreover, the efficient collisions between water drops in the 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 the 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 a larger number of aerosol particles.

The work presented here describes an extension of the libcloudph++ library
that allows us to represent the 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, if needed, further development to cover more chemical compounds and reactions.
The particle-based microphysics and chemistry scheme is used in 0-dimensional and 2-dimensional modelling set-ups.
The former set-up tests the new scheme against the
previous numerical studies that used moving-bin microphysics and aqueous-phase chemistry schemes.
The latter set-up focuses on the cloud effects on the aerosol particle size distribution (cloud-aerosol processing).
Additionally, the changes in the programming interface due to the
aqueous chemistry extension are described in Sect. A in the Appendix.
Section B 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 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.

The libcloudph++ library along with the
aqueous-phase chemistry extension,
the parcel model, and the 2-D kinematic model
are released under GNU General Public License v3.0.
The version of libcloudph++ accompanying this publication
is tagged as “2.0.0” at the project repository
and is also available as an electronic Supplement to this paper.
libcloudph++ and the 2-D slice model are
available at https://github.com/igfuw/libcloudphxx
(Arabas et al., 2015)
and the parcel model is available at
https://github.com/igfuw/parcel (last access: 28 August 2018).
The supported platforms are the following. Linux with
GNU g++, Linux with LLVM clang++, and Apple OSX with Apple clang++.
The code requires C++14 support. The compilation is
tested using the Travis continuous integration framework.

The programming interface of the particle-based microphysics scheme of libcloudph++ is
presented in Sect. 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 the C++ programming language, see Stroustrup (2013).

The aqueous chemistry module
is implemented as an optional extension to the particle-based microphysics scheme
in libcloudph++.
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 the structure named lgrngn::opts_t.
The 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 Sect. 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 the lgrngn::opts_init_t structure.
The chemistry module adds three fields to this structure:
(i) the Boolean chem_switch field that enables memory allocation for additional
variables needed for chemistry representation,
(ii) the integer sstp_chem field that defines the number of sub-steps
to be carried out in aqueous chemistry calculations, and
(iii) the real_tchem_rho field that defines 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 the 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 the code listing in Fig. A4.
The second template parameter of that structure
selects between CPU and GPU calculations
(Arabas et al., 2015).
The initialization, time stepping, and output from the particle-based scheme are done
using the methods of the lgrngn::particles_t structure.
Their signatures are provided in the code listing in
Fig. A4.

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 vapour 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 Sect. 4 none are necessary,
whereas for the 2-dimensional kinematic model from Sect. 5
two arguments are specified in order to describe the velocity field.
The last argument of init() is a map with keys from the 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.

During time stepping, the particle-based scheme computations are performed by the 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 the step_sync() method is extended
by the chemistry module.
Similar to the init() method,
a map linking the 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 the step_async method is not changed by the new
chemistry module.

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 the diag_chem method to them
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 Sect. 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.

AJ developed the aqueous-phase
chemistry extension of the libcloudph++ and carried
out the computations.
Both AJ and HP were involved in the discussion of the
results and in the
process of writing and editing the manuscript.

The work was funded by
Poland's National Science Centre (Narodowe Centrum Nauki),
grant agreements nos. 2012/06/M/ST10/00434 and 2014/15/N/ST10/05143.
We would like to thank Sonia Kreidenweis for her help
when designing and implementing the aqueous chemistry scheme
and Shin-ichiro Shima for suggesting studying the Hoppel gap formation
using the super-droplet method.
We would like to thank the two anonymous referees and Sylwester Arabas
for their comments that greatly improved the manuscript.
We would also like to thank
Travis CI and GitHub for providing their platforms free of charge
for open-source projects.

Hoffmann, F.: On the limits of Köhler activation theory: how do collision
and coalescence affect the activation of aerosols?, Atmos. Chem. Phys., 17,
8343–8356, https://doi.org/10.5194/acp-17-8343-2017, 2017. a

Riechelmann, T., Noh, Y., and Raasch, S.: A new method for large-eddy
simulations of clouds with Lagrangian droplets including the effects of
turbulent collision, New J. Phys., 14, 065008,
https://doi.org/10.1088/1367-2630/14/6/065008, 2012. a

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:
KH2O=[H+][OH-]

libcloudph++ is a free and open-source library of schemes representing cloud microphysics (e.g. condensation of water vapour into cloud droplets, collisions between water drops, precipitation) in numerical models. This work adds new schemes that represent aqueous chemical reactions in water drops. The schemes focus on the oxidation of SO2 by O3 and H2O2. The libcloudph++ is now capable of resolving the changes in aerosol sizes caused by both collisions between water drops and aqueous oxidation.

libcloudph++ is a free and open-source library of schemes representing cloud microphysics (e.g....