Retrieval of process rate parameters in the general dynamic equation for aerosols using Bayesian state estimation: BAYROSOL1.0

The uncertainty in the radiative forcing caused by aerosols and its effect on climate change calls for research to improve knowledge of the aerosol particle formation and growth processes. While experimental research has provided a large amount of high-quality data on aerosols over the last 2 decades, the inference of the process rates is still inadequate, mainly due to limitations in the analysis of data. This paper focuses on developing computational methods to infer aerosol process rates from size distribution measurements. In the proposed approach, the temporal evolution of aerosol size distributions is modeled with the general dynamic equation (GDE) equipped with stochastic terms that account for the uncertainties of the process rates. The time-dependent particle size distribution and the rates of the underlying formation and growth processes are reconstructed based on time series of particle analyzer data using Bayesian state estimation – which not only provides (point) estimates for the process rates but also enables quantification of their uncertainties. The feasibility of the proposed computational framework is demonstrated by a set of numerical simulation studies.

1 Introduction

Aerosols scatter and absorb solar radiation and affect the permeability of the atmosphere to solar energy (the direct effect). In addition, aerosol particles act as seeds for cloud droplets (cloud condensation nuclei, CCN) and, thus, influence the properties of clouds (the indirect effect; e.g., Ramanathan et al., 2001). Worldwide, particulate air pollutants are also responsible for up to 7 million premature deaths per year (WHO 2014). The Intergovernmental Panel on Climate Change (IPCC Assessment Report 2013 by Stocker et al., 2014, and IPCC Assessment Report 2014 by Pachauri et al., 2014) recognized the uncertainty in the radiative forcing caused by aerosols as the main individual factor limiting the scientific understanding of future and past climate changes.

Some of the uncertainty is caused by the fact that the initial stages of the new particle formation (NPF) processes in the atmosphere are still not completely known. It has been known for close to 100 years that photochemically driven NPF may occur in the atmosphere (Aitken, 1889). The fact that it occurs regularly throughout the troposphere has, however, only become clear during the last 15 years or so (Kulmala et al., 2004). Partly because of this, the systematic development of parameterizations describing tropospheric NPF as an active research topic has only just begun. Studies with these models suggest that tropospheric NPF may have significant effects on CCN and, in turn, on the general global cloud albedo effects of atmospheric aerosols (Merikanto et al., 2010). One challenge in estimating the anthropogenic aerosol effect on climate is the need to know the preindustrial conditions and dynamics, which is the baseline of forcing estimations (Hansen et al., 1981). Very recently, Gordon et al. (2016) made new estimations of anthropogenic aerosol radiative forcing by assuming pure biogenic particle formation, as suggested by Kirkby et al. (2016) based on the Cosmics Leaving Outdoor Droplets (CLOUD) experiment at CERN. The increased particle formation rates under preindustrial conditions increased the CCN conditions as well as the cloud albedo – resulting in a 0.22 W m⁻² decrease in anthropogenic aerosol forcing. As the NPF treatment in global climate models is typically based on parameterizing particle formation and growth rates (GRs) from chamber experiments such as CLOUD or from field measurements, it is of great importance to analyze the data with care and pay attention to uncertainties.

Most of the research dealing with NPF event analysis has used similar methodology as outlined in the “protocol” described by Kulmala et al. (2012) and Dada et al. (2020). Particle formation rates, at the detection limit of the instrument (typically in the range from 2 to 3 nm), are usually estimated from the time evolution of either the total number concentration or the number concentration below a certain size (e.g., 20 nm), correcting for coagulation loss and condensation growth with very simple and approximate balance equations. For condensational growth rate, there have been three main approaches: (1) fitting the growing nucleation mode with a lognormal function, and the growth of the mode GR is defined as the growth of the geometric mean size of the mode (Leppä et al., 2011); (2) the so-called “maximum concentration method” in which GR is estimated from the times when each measurement channel of the instrument reaches its maximum concentration (Lehtinen and Kulmala, 2003); and (3) the so-called “appearance time method” in which GR is estimated from concentration rise times of each channel (Lehtipalo et al., 2014). Methods 1 and 2 are applicable to cases in which there is a clear nucleation mode growing, as is the case, for example, for the multitude of events analyzed from the Hyytiälä forest station in Finland (Maso et al., 2005). For chamber experiments, where the aerosol size distribution approaches steady state (Dada et al., 2020), these two approaches cannot be used. Method 3 can then be applied to the transition stage of the dynamics, before steady-state is reached. All of these methods suffer from disturbances by other aerosol microphysical processes (e.g., coagulation and deposition) and, in addition, cannot be used to estimate uncertainties related to GR. Deposition rates are typically estimated by targeted experiments with either several different experiments with different monodisperse aerosol or in the absence of vapors and with low enough concentrations that other microphysical processes do not affect the estimation.

The last decade has been a huge leap forward in atmospheric NPF research. Instrument development, especially advances in particle detection efficiency and mass spectrometry, have allowed us to measure details of the dynamics of the smallest clusters (e.g., Almeida et al., 2013). At the same time, however, potentially superior advanced data analysis methods have not been used. Instead, NPF and particle growth rates have been analyzed with the abovementioned very simple regression or balance equation approaches (an overview of the typical methods is given in Kulmala et al., 2012), suffering from potentially crude approximations and permitting no proper estimation of the uncertainties. It is likely that there are significant inaccuracies in quantities such as particle formation and growth rates estimated previously, and at least their uncertainties have typically not been quantified. Very recent studies by Kürten et al. (2018) have already shown that the difference between nucleation rates estimated by fitting a sophisticated aerosol model to data and a “traditional” simple method can be as large as a factor of 10.

Studies on applying computational inversion methods to estimating the most important quantities of interest with respect to particle fate and effects in the atmosphere, the formation and growth rates, are rare. Lehtinen et al. (2004) applied simple least-squares-based optimization of aerosol microphysics to measured data. This method was later improved (more processes, less assumptions) by Verheggen and Mozurkewich (2006) and Kuang et al. (2012). These studies, however, did not address the uncertainties in the estimated parameters. Henze et al. (2004) used the method of adjoint equations to estimate condensation rates based on measured evolution of the aerosol size distribution. Sandu et al. (2005) presented adjoint equations of the complete aerosol general dynamic equation (GDE) which could be a basis for data-assimilation of aerosol dynamics. We are, however, not aware of the methodology being used later.

In the statistical (Bayesian) framework of inverse problems (Kaipio and Somersalo, 2006), the uncertainties of the model quantities are modeled statistically, and it offers an approach to uncertainty quantification, in addition to parameter estimation. In a time-invariant case, the Bayesian approach was adopted to estimate aerosol size distributions by Voutilainen et al. (2001). Thus far, the only works where the statistical approach has been taken to inverse problems in aerosol size distribution dynamics are those on parameter estimation in aggregation–fragmentation models (Ramachandran and Barton, 2010; Bortz et al., 2015; Shcherbacheva et al., 2020), estimating the size distribution evolution using Kalman filtering (Voutilainen and Kaipio, 2002; Viskari et al., 2012) and estimating evaporation rates using a Markov chain Monte Carlo method (Kupiainen-Määttä, 2016). However, the statistical inversion framework (and Bayesian state estimation in particular) has been applied to several other problems which are mathematically similar to parameter estimation in aerosol dynamics. Unknown coefficients have been estimated in, for example, Fokker–Planck equations (Banks et al., 1993; Dimitriu, 2002), age-structured population dynamics models (Rundell, 1993; Cho and Kwon, 1997) as well as algal and phytoplankton aggregation models (Ackleh, 1997; Ackleh and Miller, 2018).

In this paper, we approach the problem of estimating unknown rate parameters in the aerosol general dynamic equation in the framework of Bayesian state estimation. We model the discretized particle size distribution as well as the unknown nucleation, growth and deposition rates in GDE as multivariate random processes, and estimate them from sequential particle counter measurements using the extended Kalman filter (EKF) and fixed interval Kalman smoother (FIKS). The feasibility of these estimators to quantify the process rate parameters and their uncertainties is tested with series of numerical simulation studies.

2 Estimation of parameters in GDE

The temporal evolution of the aerosol size distribution $n=n(d_{\mathrm{p}},t)$ can be described by a population balance equation referred to as the GDE (Zhang et al., 1999; Prakash et al., 2003; Lehtinen and Zachariah, 2001; Smoluchowski, 1916; Friedlander and Wang, 1966). We write the continuous form of the GDE as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial n}{\partial t}}(d_{\mathrm{p}},t)=-\underset{\text{growth by condensation}}{\underbrace{{\displaystyle \frac{\partial g(d_{\mathrm{p}},t)n(d_{\mathrm{p}},t)}{\partial d_{\mathrm{p}}}}}}\\ {\displaystyle }& {\displaystyle }-\underset{\text{coagulation sink}}{\underbrace{n(d_{\mathrm{p}},t)\underset{d_{\mathrm{p}}^{*}}{\overset{\mathrm{\infty }}{\int }}\mathit{\beta }(d_{\mathrm{p}},s)n(s,t)\mathrm{d}s}}\\ {\displaystyle }& {\displaystyle }+\underset{\text{coagulation source}}{\underbrace{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{d_{\mathrm{p}}^{*}}{\overset{d_{\mathrm{p}}}{\int }}\mathit{\beta }\left(\sqrt[\mathrm{3}]{d_{\mathrm{p}}^{\mathrm{3}}-q^{\mathrm{3}}},q\right)n\left(\sqrt[\mathrm{3}]{d_{\mathrm{p}}^{\mathrm{3}}-q^{\mathrm{3}}},t\right)n(q,t)\mathrm{d}q}}\\ \text{(1)}& {\displaystyle }& {\displaystyle }-\underset{\text{loss by deposition}}{\underbrace{\mathit{\lambda }(d_{\mathrm{p}},t)n(d_{\mathrm{p}},t)}},\end{array}$$

where d_p is the particle diameter and t is time. Further, $g=g(d_{\mathrm{p}},t)$ denotes the condensational growth rate, $\mathit{\beta }=\mathit{\beta }(s,d_{\mathrm{p}}-s)$ is the coagulation frequency and $\mathit{\lambda }=\mathit{\lambda }(d_{\mathrm{p}},t)$ is the deposition rate.

The boundary conditions consist of fluxes of particles in and out of the considered size range $[d_{\mathrm{p}}^{\text{min}},d_{\mathrm{p}}^{\mathrm{\infty }}]$. In reality, the formation of particles occurs at very low size (typically at the 1.5–2 nm scale) by nucleation; the theoretical size at which molecule clusters start being stable is referred to as the critical size, and it is denoted $d_{\mathrm{p}}^{*}$ in Eq. (1). Note that we do not include an explicit source term in Eq. (1) to model the true nucleation (spawning particles out of vapor) as we are considering typical size ranges for particle mobility measurements (differential mobility particle sizer, DMPS, or scanning mobility particle sizer, SMPS) which are above the nucleation size. In practice, this means that $d_{\mathrm{p}}^{\text{min}}$ is always chosen to be sufficiently larger than the critical size $d_{\mathrm{p}}^{*}$. The appearance of new particles to the measurement range then occurs by condensational growth of freshly nucleated particles from below the measurement range. This process, sometimes called apparent particle formation (e.g., Lehtinen et al., 2007), is conveniently treated as a particle-concentration flux, in particle-size space (${\mathrm{cm}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$), as a boundary condition for the GDE. Hence, the nucleation J=J(t) is identified with the flux of particles to the smallest size class by condensation – that is,

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }g\left(d_{\mathrm{p}}^{\text{min}},t\right)n\left(d_{\mathrm{p}}^{\text{min}},t\right)=J\left(t\right).\end{array}$$

Similarly, for the outward flux of particles at the largest size class, we write

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }g\left(d_{\mathrm{p}}^{\mathrm{\infty }},t\right)n\left(d_{\mathrm{p}}^{\mathrm{\infty }},t\right)=\mathrm{0},\end{array}$$

which states that the growth of particles to a size exceeding $d_{\mathrm{p}}^{\mathrm{\infty }}$, the upper limit of the size range, is negligible.

The time- and/or size-dependent parameters g, β, λ and J characterize the microphysical properties of aerosols: these parameters, along with boundary conditions of the GDE, completely determine the evolution of the aerosol size distribution. However, the process rate parameters are usually not known. In this paper, we aim at estimating the growth, deposition and nucleation rates (g, λ and J, respectively) based on particle size distribution measurements. For the coagulation coefficient β, a fixed (known/approximate) value will be used.

The analysis proposed in this paper is applicable to both differential mobility particle sizer (DMPS) and scanning mobility particle sizer (SMPS) measurements. For the rest of the paper, however, we refer to measurement modality as SMPS because the number of particle size classes is relatively high in the numerical example cases.

An SMPS measurement is vector y^k∈ℝ^M that represents an indirect observation of the particle size density n(d_p,t_k), corrupted by Poisson distributed noise:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{\mathit{y}}^k={\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{y}}}^k}{V}},\text{ s.t. }{\stackrel{\mathrm{̃}}{\mathit{y}}}^k\sim \text{Poisson}\left(V{\mathit{z}}^k\right),\text{ and }{\mathit{z}}^k=\mathcal{H}n\left(d_{\mathrm{p}},t_k\right),\end{array}$$

where V is a constant (the effective volume of the sample in the condensation particle counter), ℋ is a device-dependent linear operator and M is the number of channels in the particle counter. In the following, we assume that the rate of time evolution is negligible compared with the time required to measure M channels, or one frame, with an SMPS. We denote the measurement time of the kth frame by t_k.

As SMPS measurements explicitly depend on the size density n alone and not on g, λ and J, a measurement y^k corresponding to a single time instant does not carry enough information to estimate these parameters. However, as g, λ and J determine the temporal evolution of n, it might be possible to estimate them on the basis of a sequence of measurements y^k corresponding to a set of time instants t_k, $k=\mathrm{1},\mathrm{\dots },K$.

In this section, we formulate the problem of estimating the time-dependent size density n and the process rate parameters g, λ and J as a Bayesian state estimation problem. To this end, we first discretize the GDE with respect to size and time, and write it in a stochastic form in order to model its uncertainties. We also model the process rate parameters as discrete-time stochastic processes. This formulation allows us to express the following questions in the Bayesian framework (Gelb, 1974):

• What are the expected values of n(d_p,t_k), g(d_p,t_k), λ(d_p,t_k) and J(t_k) at each time t_k given a set of measurements ${\mathbf{Y}}^{\mathrm{\ell }}=\mathit{\{}{\mathit{y}}^{\mathrm{1}},\mathrm{\dots },{\mathit{y}}^{\mathrm{\ell }}\mathit{\}}$ corresponding to discrete times $t_{\mathrm{1}},\mathrm{\dots },t_{\mathrm{\ell }}$?

• How large are the uncertainties of the estimated quantities?

The state estimation problems are referred to as prediction, filtering and smoothing depending on whether ℓ<k, ℓ=k or ℓ>k, respectively. While filtering is a suitable choice for online monitoring and control problems, smoothing is usually a preferable choice when estimates are not needed online; the smoother estimates also utilize the future observations ${\mathit{y}}^{k+\mathrm{1}},\mathrm{\dots },{\mathit{y}}^{\mathrm{\ell }}$ for the estimate corresponding to time t_k.

The latter question refers to posterior uncertainties – that is, uncertainties of the quantities given the measurements Y^ℓ. In the Bayesian framework, these uncertainties can be quantified by computing, for example, posterior variances and credible intervals of the parameters.

In the simplest special case, where the evolution model and observation model are linear with respect to all parameters and all error terms are additive and Gaussian, the Bayesian filtering and smoothing problems can be solved by Kalman filter and Kalman smoother recursions, respectively (Kalman, 1960; Gelb, 1974). In general cases, where models are nonlinear or non-Gaussian, only approximate solutions are available. In principle, the best approximations of the posterior estimates are obtained with sequential Monte Carlo methods, known as particle filters/smoothers (Särkkä, 2013). However, these methods are limited to small dimensional cases, because of the high computational burden. Computationally more efficient approximations include ensemble, unscented and extended Kalman filters/smoothers (Särkkä, 2013). In this paper, we choose the extended Kalman filter (EKF) and fixed interval Kalman smoother (FIKS), but we note that the other filters and smoothers developed for nonlinear state estimation are applicable as well. In the next section, the feasibility of the EKF and FIKS for the GDE parameter estimation problem will be tested numerically.

2.1 Evolution model

2.1.1 Discretized, stochastic GDE

To approximate the GDE (1) numerically, we partition the particle size variable d_p into Q intervals (or bins) Ω_i of widths Δd_i, $i=\mathrm{1},\mathrm{\dots },Q$. The discrete instants in time in the temporal discretization are denoted by t_k, $k=\mathrm{1},\mathrm{\dots },T$, and the differences between consecutive times are denoted by $\mathrm{\Delta }{t}^k=t_{k+\mathrm{1}}-t_k$. We denote the number concentration of particles corresponding to the ith bin at time t_k by $N_i^k$ – that is, $N_i^k={\int }_{{\mathrm{\Omega }}_i}n(d_{\mathrm{p}},t_k)\mathrm{d}{d}_{\mathrm{p}}$, and a vector consisting of particle concentrations in all Q bins at time t_k by N^k, that is, ${\mathit{N}}^k=[N_{\mathrm{1}}^k,\mathrm{\dots },N_Q^k{]}^{\mathrm{T}}$. We discretize the condensation and deposition rates accordingly, and write ${\mathit{g}}^k=[g_{\mathrm{1}}^k,\mathrm{\dots },g_Q^k{]}^{\mathrm{T}}$, ${\mathit{\lambda }}^k=[{\mathit{\lambda }}_{\mathrm{1}}^k,\mathrm{\dots },{\mathit{\lambda }}_Q^k{]}^{\mathrm{T}}$. Further, the nucleation J is discretized with respect to time: J^k denotes the nucleation rate at time t_k.

Using Euler's method for time integration and first-order upwind differencing for the condensation terms, we get the discrete-time evolution model for the particle number concentrations in bins Ω_i (Korhonen et al., 2004):

$$\begin{array}{ll}\text{(5)}& {\displaystyle }& {\displaystyle }{N}_{\mathrm{1}}^{k+\mathrm{1}}=N_{\mathrm{1}}^k+\mathrm{\Delta }{t}^k\left(J^k-\left({\displaystyle \frac{g_{\mathrm{1}}^k}{\mathrm{\Delta }{d}_{\mathrm{1}}}}+{\mathit{\lambda }}_{\mathrm{1}}\right)N_{\mathrm{1}}^k\right)+C_{\mathrm{1}}\left({\mathit{N}}^k\right){\displaystyle }& {\displaystyle }{N}_i^{k+\mathrm{1}}=N_i^k+\mathrm{\Delta }{t}^k\left({\displaystyle \frac{g_{i-\mathrm{1}}^k}{\mathrm{\Delta }{d}_{i-\mathrm{1}}}}{N}_{i-\mathrm{1}}^k-\left({\displaystyle \frac{g_i^k}{\mathrm{\Delta }{d}_i}}+{\mathit{\lambda }}_i\right)N_i^k\right)\\ \text{(6)}& {\displaystyle }& {\displaystyle }+C_i\left({\mathit{N}}^k\right),{\displaystyle }& {\displaystyle }\text{for all }\mathrm{1}<i\le Q,\mathrm{1}\le k\le T,\end{array}$$

where C_i(N^k) is a nonlinear coagulation term; for details, see Lehtinen and Zachariah (2001). Here, the choice for approximating the derivative in the growth term in the discretization of the GDE is made, assuming that the particle growth rate is positive. The modification to cases where the aerosols evaporate, i.e., where growth rate is negative, is straightforward. Such a time discretization scheme can become unstable; however, it is possible to apply a different method, e.g., implicit Euler or Crank–Nicholson (Trangenstein, 2013). Equations (5) and (6) can be written equivalently in vector form as

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{\mathit{N}}^{k+\mathrm{1}}=\mathbf{A}\left({\mathit{g}}^k,{\mathit{\lambda }}^k\right){\mathit{N}}^k+s\left(J^k\right)+C\left(\mathit{\beta },{\mathit{N}}^k\right).\end{array}$$

Here, $\mathbf{A}=A({\mathit{g}}^k,{\mathit{\lambda }}^k)$ is a sparse matrix consisting of elements A_ij:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{A}_{ij}({\mathit{g}}^k,{\mathit{\lambda }}^k)=\left\{\begin{array}{ll}\mathrm{1}-\mathrm{\Delta }{t}^k\left(\frac{g_i^k}{\mathrm{\Delta }{d}_i}+{\mathit{\lambda }}_i\right)& i=j\\ \mathrm{\Delta }{t}^k\frac{g_{i-\mathrm{1}}^k}{\mathrm{\Delta }{d}_{i-\mathrm{1}}}& i=j+\mathrm{1}>\mathrm{1}\\ \mathrm{0}& \text{otherwise},\end{array}\right.\end{array}$$

$\mathit{s}=s\left(J^k\right)\in {\mathbb{R}}^Q$ is a vector of the form $s\left(J^k\right)=[\mathrm{\Delta }{t}^k{J}^k,\mathrm{0},\mathrm{\dots },\mathrm{0}{]}^{\mathrm{T}}$, and the nonlinear term C(β,N^k) is defined accordingly.

Finally, we complement the discretized GDE with a stochastic term ϵ_k∈ℝ^Q to account for modeling errors caused by factors such as discretization and uncertainties of the boundary conditions, and we write

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{\mathit{N}}^{k+\mathrm{1}}=f\left({\mathit{N}}^k\right)+{\mathit{ϵ}}_k,\end{array}$$

where $f:{\mathbb{R}}^Q\to {\mathbb{R}}^Q$ is of the form $f\left({\mathit{N}}^k\right)=\mathbf{A}({\mathit{g}}^k,{\mathit{\lambda }}^k){\mathit{N}}^k+s\left(J^k\right)+C(\mathit{\beta },{\mathit{N}}^k)$. In this paper, the stochastic state noise term ϵ_k is modeled as Gaussian, ${\mathit{ϵ}}_k\sim \mathcal{N}(\mathit{0},{\mathbf{\Gamma }}_{\mathit{ϵ}})$, where Γ_ϵ is the covariance matrix of ϵ_k. Equation (9) forms a discretized, stochastic evolution model for the particle number concentration N.

2.1.2 Models for the parameters of the GDE

All the unknown parameters of the GDE (g, λ and J) are known to be nonnegative. For this reason, we reparametrize these quantities by writing $g_i^k=P_g\left({\mathit{\xi }}_{g,i}^k\right)$, ${\mathit{\lambda }}_i^k=P_{\mathit{\lambda }}\left({\mathit{\xi }}_{\mathit{\lambda },i}^k\right)$ and $J^k=P_J\left({\mathit{\xi }}_J^k\right)$, where ξ_g,i, ξ_λ,i and ξ_J are the (unconstrained) parameters, and mappings P_g, P_λ, P_J have the form of a so-called softplus function (Dugas et al., 2001) $P_{\mathit{\phi }}:\mathbb{R}\to {\mathbb{R}}^{+}$:

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{P}_{\mathit{\phi }}\left({\mathit{\xi }}_{\mathit{\phi }}^k\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}\log \left(\mathrm{1}+e^{\mathit{\alpha }{\mathit{\xi }}_{\mathit{\phi }}^k}\right).\end{array}$$

Note that this parameterization only constrains the parameter from below, but if an upper boundary is known, the function in Eq. (10) could be changed in favor of, for example, a logistic function which allows for both constraints from below and from above. We denote the vectors consisting of all condensation and deposition parameters at time t_k by ${\mathit{\xi }}_g^k=[{\mathit{\xi }}_{g,\mathrm{1}}^k,\mathrm{\dots },{\mathit{\xi }}_{g,Q}^k{]}^{\mathrm{T}}$ and ${\mathit{\xi }}_{\mathit{\lambda }}^k=[{\mathit{\xi }}_{\mathit{\lambda },\mathrm{1}}^k,\mathrm{\dots },{\mathit{\xi }}_{\mathit{\lambda },Q}^k{]}^{\mathrm{T}}$, respectively.

As the process rates g, λ and J in the GDE are time-varying, the state estimation requires one to model their time dependence. In this paper, we model ξ_g,i, ξ_λ,i and ξ_J as either first-order Markov processes,

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }{\mathit{\xi }}_{\mathit{\phi }}^{k+\mathrm{1}}={\mathbf{\Psi }}_{\mathit{\phi }}{\mathit{\xi }}_{\mathit{\phi }}^k+{\mathit{\eta }}_{\mathit{\phi }}^k,\end{array}$$

where Ψ_φ is a diagonal matrix Ψ_φ=r_φI, and $\mathrm{0}<r_{\mathit{\phi }}<\mathrm{1}$, or as second-order Markov processes,

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{\mathit{\xi }}_{\mathit{\phi }}^{k+\mathrm{2}}={\mathbf{\Psi }}_{\mathit{\phi }}^{\mathrm{1}}{\mathit{\xi }}_{\mathit{\phi }}^{k+\mathrm{1}}+{\mathbf{\Psi }}_{\mathit{\phi }}^{\mathrm{2}}{\mathit{\xi }}_{\mathit{\phi }}^k+{\mathit{\eta }}_{\mathit{\phi }}^k,\end{array}$$

with ${\mathbf{\Psi }}_{\mathit{\phi }}^{\mathrm{1}}=r_{\mathit{\phi }}^{\mathrm{1}}\mathbf{I}$ and ${\mathbf{\Psi }}_{\mathit{\phi }}^{\mathrm{2}}=r_{\mathit{\phi }}^{\mathrm{2}}\mathbf{I}$. In both models, ${\mathit{\eta }}_{\mathit{\phi }}^k$ stands for Gaussian noise ${\mathit{\eta }}_{\mathit{\phi }}^k\sim \mathcal{N}(\mathbf{0},{\mathbf{\Gamma }}_{{\mathit{\eta }}_{\mathit{\phi }}^k})$. To simplify the following description, we assume all the models to be of the form shown in Eq. (11), but we note that the extension to second-order models is straightforward: higher-order Markov models can be converted into the form of first-order Markov models by augmenting the state variables corresponding to more than one time instant into a single vector. The second-order models are suitable for some of the quantities in GDE, because they imply temporal smoothness of those processes (Kaipio and Somersalo, 2006). The specific choices of the state models and their parameters are discussed in Sect. 3 and specified in Appendix B.

The evolution models, such as Eqs. (11) and (12) can be argued to be unrealistic, as they are not based on physics. The understanding, however, is that if the (co)variances of the driving noise processes ${\mathit{\eta }}_{\mathit{\phi }}^k$ are set high enough, such models are feasible in the sense that the actual ${\mathit{\xi }}_{\mathit{\phi }}^{k+\mathrm{1}}-{\mathbf{\Psi }}_{\mathit{\phi }}{\mathit{\xi }}_{\mathit{\phi }}^k$ are well supported by the modeled distribution of ${\mathit{\eta }}_{\mathit{\phi }}^k$. Systematic (state-space identification) approaches exist that allow one to test the model's feasibility against the driving noise distribution (variances) (Gelb, 1974). In the next section, we test the state estimation based on the above models in cases where the true evolution of the quantities is not of the form of Markov models.

2.1.3 Augmented evolution model for n, g, λ and J

To complete the evolution model, we define an augmented state variable X^k,

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{\mathit{X}}^k=\left[\begin{array}{l}{\mathit{N}}^k\\ {\mathit{\xi }}_g^k\\ {\mathit{\xi }}_{\mathit{\lambda }}^k\\ {\mathit{\xi }}_J^k\end{array}\right],\end{array}$$

and combine the evolution models written in Sects. 2.1.1 and 2.1.2, yielding

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left[\begin{array}{l}{\mathit{N}}^{k+\mathrm{1}}\\ {\mathit{\xi }}_g^{k+\mathrm{1}}\\ {\mathit{\xi }}_{\mathit{\lambda }}^{k+\mathrm{1}}\\ {\mathit{\xi }}_J^{k+\mathrm{1}}\end{array}\right]=\\ {\displaystyle }& {\displaystyle }\left[\begin{array}{cccc}A\left(P_g\right({\mathit{\xi }}_{g,i}^k),P_{\mathit{\lambda }}({\mathit{\xi }}_{\mathit{\lambda },i}^k\left)\right)& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{\Psi }}_g& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{\Psi }}_{\mathit{\lambda }}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathbf{\Psi }}_J\end{array}\right]\times \left[\begin{array}{l}{\mathit{N}}^k\\ {\mathit{\xi }}_g^k\\ {\mathit{\xi }}_{\mathit{\lambda }}^k\\ {\mathit{\xi }}_J^k\end{array}\right]\\ \text{(14)}& {\displaystyle }& {\displaystyle }+\left[\begin{array}{c}s\left(P_J\right({\mathit{\xi }}_J^k\left)\right)\\ \mathit{0}\\ \mathit{0}\\ \mathit{0}\end{array}\right]+\left[\begin{array}{c}C(\mathit{\beta },{\mathit{N}}^k)\\ \mathit{0}\\ \mathit{0}\\ \mathit{0}\end{array}\right]+\left[\begin{array}{c}{\mathit{ϵ}}^k\\ {\mathit{\eta }}_g^k\\ {\mathit{\eta }}_{\mathit{\lambda }}^k\\ {\mathit{\eta }}_J^k\end{array}\right]\end{array}$$

or

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }{\mathit{X}}^{k+\mathrm{1}}=F\left({\mathit{X}}^k\right)+{\mathit{w}}^k.\end{array}$$

This is the evolution model for the augmented state variable X^k, which includes not only the number concentrations of the size sections but also the unknown process rates. Next, we write the observation model in terms of X^k, and then, in Sect. 2.3, we apply Bayesian state estimation to infer X^k, $k=\mathrm{1},\mathrm{\dots },K$ based on sequential SMPS measurements.

2.2 Observation model

An SMPS consists of a differential mobility analyzer (DMA), which classifies charged particles based on their mobility in an electric field, and a condensation particle counter (CPC) where the classified particles are grown to sizes detectable, for example, optically. All particle counters provide only discrete, indirect and noisy data on the particle size distributions. Mathematically, each channel in a particle counter gives data that corresponds to a convolution/projection of the particle size distribution onto a space spanned by device-specific kernel functions; moreover, the counter data are corrupted by Poisson distributed noise.

2.2.1 SMPS transfer function

The output of each DMA is the detected number concentration in a size class; in channel i corresponding to a discrete time index k, it is of the form

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }{z}_i^k={\displaystyle \frac{\mathrm{1}}{V}}\underset{t_{\mathrm{0}}+(k-\mathrm{1})\mathrm{\Delta }t}{\overset{t_{\mathrm{0}}+k\mathrm{\Delta }t}{\int }}{\mathit{\varphi }}_a\left(t\right)\underset{{\mathit{\omega }}_i}{\int }{\mathrm{\Psi }}_i\left(d_{\mathrm{p}}\right)n(d_{\mathrm{p}},t)\mathrm{d}{d}_{\mathrm{p}}\mathrm{d}t,\end{array}$$

where Ψ_i(d_p) is a time-invariant kernel function; ω_i is the support of Ψ_i(d_p) – that is, the set of points where Ψ_i(d_p) is nonzero; t₀ is the initial time; and Δt is the duration of counting particles in the CPC for a single channel. Further, V is the volume of the aerosol sample that passes through the CPC counter with a detector–sample flow rate ϕ_a(t) in the period of time Δt – that is, $V={\int }_{t_{\mathrm{0}}+(k-\mathrm{1})\mathrm{\Delta }t}^{t_{\mathrm{0}}+k\mathrm{\Delta }t}{\mathit{\varphi }}_a\left(t\right)\mathrm{d}t$.

2.2.2 CPC counting model

The measurement data of the ith channel of an SMPS, $y_i^k$, consist of Poisson-distributed counts given by the CPC:

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }{y}_i^k={\displaystyle \frac{{\stackrel{\mathrm{̃}}{y}}_i^k}{V}},\text{with }{\stackrel{\mathrm{̃}}{y}}_i^k\sim \text{Poisson}\left(V{z}_i^k\right),\end{array}$$

where $V\simeq \mathrm{\Delta }t{\mathit{\varphi }}_a(t_{\mathrm{0}}+k\mathrm{\Delta }t)$ is the volume of sample used in the CPC for counting. In the numerical studies of this paper, we use this model when simulating the measurement data. In state estimation, however, we approximate the Poisson distributed observations as Gaussian:

$$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle }{y}_i^k\sim \mathcal{N}\left(z_i^k,{\mathit{\gamma }}_i\right),\end{array}$$

where γ_i is the approximate variance of the noise. In this paper, we use the same approximation as in Voutilainen et al. (2000) and write ${\mathit{\gamma }}_i^k=\frac{y_i^k}{V}$.

By discretizing the SMPS model (Eq. 16) and using the above Gaussian approximation of the noise, we write an observation model of the form

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }{\mathit{y}}^k=\overline{\mathbf{H}}{\mathit{N}}^k+{\stackrel{\mathrm{̃}}{\mathit{e}}}^k+{\mathit{\iota }}^k,\end{array}$$

where ${\mathit{y}}^k=[y_{\mathrm{1}}^k,\mathrm{\dots },y_M^k{]}^{\mathrm{T}}$, $\overline{\mathbf{H}}$ is an observation matrix and ${\stackrel{\mathrm{̃}}{\mathit{e}}}^k$ is the Gaussian observation noise ${\stackrel{\mathrm{̃}}{\mathit{e}}}^k\sim \mathcal{N}(\mathit{0},{\mathbf{\Gamma }}_{\stackrel{\mathrm{̃}}{\mathit{e}}}^k)$. Here, the covariance of the observation noise is of the diagonal form ${\mathbf{\Gamma }}_e^k=\text{diag}({\mathit{\gamma }}_{\mathrm{1}}^k,\mathrm{\dots },{\mathit{\gamma }}_M^k)$. The additional noise term, ι^k is included in Eq. (19) in order to account for errors caused by the discretization of the measurement operator in Eq. (16). Here, ι^k is simply approximated as being Gaussian distributed with zero-mean ${\mathit{\iota }}^k\sim \mathcal{N}(\mathit{0},{\mathbf{\Gamma }}_{\mathit{\iota }}^k)$ whose components are calculated as mutually independent; hence, the covariance matrix ${\mathbf{\Gamma }}_{\mathit{\iota }}^k$ is of diagonal form. Furthermore, the approximation error term ι^k is assumed to be independent of the counting noise ${\stackrel{\mathrm{̃}}{\mathit{e}}}^k$; thus, the total error ${\mathit{e}}^k={\stackrel{\mathrm{̃}}{\mathit{e}}}^k+{\mathit{\iota }}^k\sim \mathcal{N}(\mathit{0},{\mathbf{\Gamma }}_e^k)$, where ${\mathbf{\Gamma }}_e^k={\mathbf{\Gamma }}_{\stackrel{\mathrm{̃}}{\mathit{e}}}^k+{\mathbf{\Gamma }}_{\mathit{\iota }}^k$. For a more rigorous approach to handling modeling errors, we refer to the book by Kaipio and Somersalo (2006).

Finally, the model (Eq. 19) can be written in terms of the state variable X^k defined in Eq. (13),

$$\begin{array}{}\text{(20)}& {\displaystyle }{\displaystyle }{\mathit{y}}^k=\mathbf{H}{\mathit{X}}^k+{\mathit{e}}^k,\end{array}$$

where H is a block matrix of the form $\mathbf{H}=[\overline{\mathbf{H}},\mathbf{0},\mathbf{0},\mathbf{0}]$.

2.3 State estimation

The nonlinear evolution model (Eq. 15) and the linear observation model (Eq. 20) form a system

$$\begin{array}{}\text{(21)}& {\displaystyle }& {\displaystyle }{\mathit{X}}^{k+\mathrm{1}}=F\left({\mathit{X}}^k\right)+{\mathit{w}}^k\text{(22)}& {\displaystyle }& {\displaystyle }{\mathit{y}}^k=\mathbf{H}{\mathit{X}}^k+{\mathit{e}}^k,\end{array}$$

referred to as the state-space representation. The system is stochastic due to the state noise and observation noise processes, w^k and e^k, respectively. In addition, we model the initial state X⁰ as a Gaussian random variable ${\mathit{X}}^{\mathrm{0}}\sim \mathcal{N}({\mathit{X}}^{\mathrm{0}|\mathrm{0}},{\mathbf{\Gamma }}^{\mathrm{0}|\mathrm{0}})$. Given this model, we can state the Bayesian filtering and smoothing problems as follows: form the conditional probability density of the random variable X^k, given the sequence of measurements ${\mathbf{Y}}^{\mathrm{\ell }}=\mathit{\{}{\mathit{y}}^{\mathrm{1}},\mathrm{\dots },{\mathit{y}}^{\mathrm{\ell }}\mathit{\}}$. In the extended Kalman filter and smoother, these probability densities π(X^k|Y^ℓ), are approximated by Gaussian densities – that is,

$$\begin{array}{}\text{(23)}& {\displaystyle }{\displaystyle }\mathit{\pi }\left({\mathit{X}}^k\right|{\mathbf{Y}}^{\mathrm{\ell }})\approx \mathcal{N}({\mathit{X}}^{k|\mathrm{\ell }},{\mathbf{\Gamma }}^{k|\mathrm{\ell }}),\end{array}$$

where X^k|ℓ and Γ^k|ℓ are approximations of the respective conditional expectation and covariance of X^k given Y^ℓ (Gelb, 1974).

Kalman filtering gives online estimates π(X^k|Y^k) based on the data set from the beginning up to the present time k – that is, ℓ=k, whereas in smoothing ℓ>k. In the FIKS, in particular, ℓ=K, where K is the index of the final time step. In other words, the FIKS estimate π(X^k|Y^K) at each time k is based on the entire data set from the beginning to the end of measurements. The EKF and FIKS estimates (approximate expectations X^k|ℓ and covariances Γ^k|ℓ) are given by the following forward and backward iterations (Gelb, 1974):

In Algorithms 1 and 2, ∂F^k denotes the Jacobian matrix of the nonlinear mapping F at point X^k|k. Note that the FIKS is based on a backward recursion, which starts from the filter estimate corresponding to final state: X^K|K, Γ^K|K.

3 Numerical simulations

In this section, the feasibility of the proposed estimation scheme is tested with numerical studies, where the aerosol particle evolution is simulated by numerical approximations of the GDE corresponding to a set of process rates, and where synthetic SMPS data are computed by numerical modeling of the DMA and generating the Poisson-distributed CPC data. Two type of events are considered: in cases 1 and 2, the evolution of the aerosol size distribution is governed by a nucleation event (NE) and the subsequent growth of the nucleation mode in a background of existing aerosols; in cases 3 and 4, new particles are formed in a continuous nucleation process and their further evolution is controlled by condensational growth and deposition, so that the size distribution approaches a steady state (SS). In all cases, the particle growth is dominated by condensation, and the loss of particles (caused by wall deposition, sedimentation, dilution, etc.) depends linearly on the size distribution function. Further, the coagulation kernel is chosen to have the form given in the book by Seinfeld and Pandis (2016). In these numerical studies, cases 1 and 2 (NE) qualitatively represent a typical particle formation event in the atmosphere (e.g., Hyytiälä Maso et al., 2005), whereas cases 3 and 4 (SS) represent particle formation and growth in a chamber experiment (e.g., CLOUD Lehtipalo et al., 2014).

3.1 Cases 1 and 2: nucleation event (NE)

3.1.1 NE: data simulation

In the numerical simulation study, the temporally evolving particle size distribution is synthetically generated using the (deterministic) discretized GDE model described by Eqs. (5) and (6) with predefined process rates g,λ, β and J. In the estimation, however, g,λ and J are, of course, not known.

In the NE case, the process rates g,λ and J are chosen to have the following properties: the condensational growth rate g is independent of the particle size but depends on time, whereas the loss rate λ depends on particle size but not on time. Further, the nucleation rate J is a time-dependent, continuous function which represents an NE in a time interval $[t_{\mathrm{0}}=\mathrm{5}\mathrm{h},t_{\mathrm{1}}=\mathrm{10}\mathrm{h}]$ and is zero in all other instants within the period of interest [0,15 h]. As noted above, the coagulation kernel is known, yet the coarse discretization also causes error in this term. The detailed forms of the process rates are shown in Table A1 (Appendix A).

When simulating the evolution of the particle distribution, the particle diameter $d_{\mathrm{p}}\in [\mathrm{13.85}\mathrm{nm},\mathrm{1000.0}\mathrm{nm}]$ is discretized into Q=2500 logarithmically distributed bins. Such a high size resolution is chosen in order to avoid numerical diffusion effects and to obtain a good approximation of the particle size density. The time step Δt^k in the explicit Euler time integration scheme is chosen based on the Courant–Friedrichs–Lewy condition (Courant et al., 1928; Dullemond and Dominik, 2005):

$$\begin{array}{}\text{(24)}& {\displaystyle }{\displaystyle }\mathrm{0}<{\mathrm{\Delta }}_t^k<{\displaystyle \frac{\mathrm{1}}{\underset{i}{max}\left\{\frac{g_i^k}{\mathrm{\Delta }{d}_i}+{\mathit{\lambda }}_i\right\}}}.\end{array}$$

The CFL criterion is applied in order to keep the time integration stable with respect to condensational growth and loss. The coagulation rates are not considered here because coagulation is actually a dampening mechanism that stabilizes time integration. The nucleation, growth and loss rates as well as the resulting particle size density evolution for the NE cases are illustrated in Fig. 1.

Figure 1 Cases 1 and 2 – a nucleation event, showing true (simulated) processes: (a) particle size density (measurement), (b) nucleation rate, (c) growth rate and (d) loss rate.

Figure 2 The size distribution transfer functions for (a) simulating the measurement data, and (b) the observation model used in state estimation. In the images, the horizontal axis corresponds to the size of the particle entering the device, and the vertical axis represents the channels of the SMPS. The colors represent the values of the efficiency with which a given particle size will be classified in a given bin in the particle sizer.

As explained in Sect. 2.2, the measurements consist of simulated counts modeling CPC combined with DMA. The model for the kernel functions Ψ_i(d_p) corresponds to an SMPS 3936 device, and it accounts for the most relevant effects (Millikan, 1923; Stolzenburg, 1989; Flagan, 1998; McMurry, 2000; Boisdron and Brock, 1970; Wiedensohler, 1988). We skip the details here, and only visualize the kernels, by plotting the size-distribution transfer function, or the observation matrix $\overline{\mathbf{H}}$, as a color map (Fig. 2a).

We simulate the synthetic CPC measurement $y_i^k$ corresponding to each CPC channel i at each time k by drawing samples from a Poisson distribution, given in Eq. (17) with mean $V{z}_i^k$. As the expectation of $y_i^k$ is $z_i^k$ and its variance is $z_i^k/V$, the signal-to-noise ratio (SNR) of CPC data increases with V. In order to investigate the effect of SNR to the process rate estimates, we generate the Poisson-distributed observations corresponding to two sample volumes: V=90 cm³ (Case 1: NE, high SNR) and V=0.9 cm³ (Case 2: NE, low SNR). In Case 1, the range of SNR was [0,6426]; in Case 2 it was [0,64.26].

3.1.2 NE: parameter estimation

In this section, we briefly describe the assumptions made when constructing the models used for computing the state estimates in the NE cases. The exact forms of the models as well as choices of parameters are listed in Appendix B.

For the evolution models used in state estimation, the particle size range [14.1 nm, 736.5 nm] is divided into 111 bins – that is, the size range is narrower and the discretization is significantly coarser than when simulating the data. The GDE-based, discrete stochastic state evolution model (Eq. 9) for the particle number N is written as described in Sect. 2.1.1. The covariance matrix Γ_ϵ of the stochastic term ϵ_k is chosen to be of diagonal form.

In this case study, we assume to know that the condensation rate is time dependent but size independent. Moreover, both the condensation and nucleation rates are assumed to be temporally smooth, and we model them as second-order stochastic processes (Eq. 12). The loss factor λ is assumed to be a smooth function of the particle size. Further, we write a first-order Markov model (Eq. 11) for its temporal variation – that is, although the true loss factor is time invariant, we do not assume to know this property in state estimation. This is done to study the stability of the estimation scheme: although the deposition loss is modeled as varying with time, the estimation should yield essentially time-invariant estimates.

The covariance of the initial state, Γ^0|0, is chosen to be diagonal; this signifies that the elements of X⁰ are mutually independent. Moreover, the variances of X⁰ are chosen to be relatively large in comparison with the variances of the state noise vector ϵ_k; this indicates a high uncertainty of the initial state. We note that the selection of the parameters in the stochastic terms is a crucial part of the state-space model. However, the state estimates are not extremely sensitive to these choices; choosing parameters that are of the right order of magnitude is usually enough, and as the stochastic models are written for physically relevant quantities, ballpark ranges of the parameters are often available a priori.

The size-distribution transfer function corresponding to the discretization of the particle size in estimation is illustrated in Fig. 2b. The approximate observation model is used in order to avoid so-called inverse crime, which refers to the use of unrealistically accurate models in the inversion of simulated data. A comparison between the true transfer function and the approximated one used in the estimation – mean over each discretization bins – is depicted in Fig. 2c for the third channel. Furthermore, as noted in Sect. 2.2, instead of using the Poisson model for the measurements, the observation noise is approximated as additive and Gaussian. This choice is made for computational convenience, as it allows for the direct applications of EKF and FIKS into the state-space system.

The extended Kalman filter and smoother estimates are computed using Algorithms 1 and 2, respectively. From the resulting state estimates X^k|ℓ, $\mathrm{\ell }=k,K$ the approximate posterior expectations of the processes are computed using the models $g_i^k=P_g\left({\mathit{\xi }}_{g,i}^k\right)$, ${\mathit{\lambda }}_i^k=P_{\mathit{\lambda }}\left({\mathit{\xi }}_{\mathit{\lambda },i}^k\right)$ and $J^k=P_J\left({\mathit{\xi }}_J^k\right)$ (see Sects. 2.1.2 and 2.1.3). We also compute approximate 68 % credible intervals of the estimates, by mapping the values $E\left({\mathit{\xi }}_{*,i}^k\right|{\mathbf{Y}}^{\mathrm{\ell }})\pm \sqrt{\text{var}\left({\mathit{\xi }}_{*,i}^k\right|{\mathbf{Y}}^{\mathrm{\ell }})}$, to the corresponding process rate spaces. Note, however, that due to the linearizations/Gaussian approximation behind EKF, these approximate intervals do not necessarily represent the ranges where the true parameter value lies; these are the ranges within which the true value most likely lies with an 0.68 probability. In the following, we refer to these approximate credible intervals as“posterior error intervals”.

3.1.3 NE: results and discussion

The results of Case 1 (NE, high SNR) are illustrated in Fig. 3. The figure shows the Kalman smoother estimates for the particle size density, and the EKF and FIKS estimates for the growth, apparent particle formation and loss rates as well as for two instantaneous particle size densities ($\frac{\mathrm{\Delta }N}{\mathrm{\Delta }{d}_{\mathrm{p}}}$ at 2 and 10 h). The loss rate estimates correspond to time 10 h. For the process rates and for the instantaneous size densities, the figure also shows the EKF- and FIKS-based posterior error intervals – representing the uncertainty of the estimates – as well as the true values of the corresponding quantities.

Figure 3 Case 1 – NE and high SNR. State estimates for the particle size density (a, e, f), nucleation rate (b), growth rate (c) and loss rate (d). The image in plot (a) depicts the approximate posterior expectation for the entire time-evolution of the particle size density given by FIKS, whereas plots (e) and (f) illustrate the EKF and FIKS estimates corresponding to times 2 and 10 h, respectively. In plots (b)–(f), the blue and orange lines represent the approximate posterior expectations for EKF and FIKS, respectively, and the areas shaded with light blue and orange are the respective posterior error intervals. The true values of the corresponding quantities are drawn with a green line.

The estimated particle size density (Fig. 3a) is in rather good correspondence with the true size density (Fig. 1a). However, the size density estimates corresponding to times 2 and 10 h (Fig. 3e and f, respectively) show that the peak values of the size density are somewhat underestimated by both EKF and FIKS. Moreover, in these instants, the true values of the size density are partly outside the approximate posterior error intervals. Yet the 68 % posterior error intervals do not necessarily need to contain the true values, and the error intervals seem to be slightly too narrow. This underestimation of the uncertainty is due to the linearizations/Gaussian approximation behind EKF and FIKS, as shown by Huttunen et al. (2018), and the rather simple approximations of the error in the evolution and the measurement models – underestimating the covariance in the models. The errors caused by model approximations become more influential with decreasing mean noise level – the remedy for such errors could be the Bayesian approximation error method (Kaipio and Somersalo, 2006); however, this is out of the scope of this paper.

For the process rates, the approximate posterior means given by both EKF and FIKS are relatively close to the true values (Fig. 3b–d). Overall, the FIKS estimates for the process rate parameters are more accurate than the EKF estimates – this is an expected result because FIKS utilizes the entire data set up to the end of the process, whereas EKF only uses data up to time t when estimating the variables at time t.

The filter and smoother also show differences in the posterior error intervals of the process rate parameters: the error intervals given by EKF are systematically wider than those given by FIKS. This is again an expected result because the use of the future data (FIKS) should reduce the uncertainty in the estimated quantities. Furthermore, in almost all instants in time, the process rate parameters (especially growth and apparent particle formation rates) are within the posterior error intervals. This is a desired result, as it indicates that these approximate credible intervals give realistic measures of the estimate uncertainties in these cases.

The loss rate estimate uncertainty depends strongly on the particle size: the posterior error intervals are wide in the lowest and highest size ranges and rather narrow elsewhere. The high uncertainty of the loss factor at the high particle size range d_p>400 nm is caused by the lack of data – in this size range, the particle density is nearly zero at all times; consequently, the SMPS data do not provide information on the loss factor. In the lowest size range, the width of the credible interval also depends on time. From time t = 9 h, when the nucleation event is practically over, the particle size density in the lowest size range is almost zero, again resulting in high uncertainty in the loss factor estimate. The loss rate estimate and posterior error interval resulting from the FIKS are virtually time invariant, even though those from the EKF show a strong time dependence.

Figure 4 Case 2 – NE and low SNR. State estimates for the particle size density (a, e, f), nucleation rate (b), growth rate (c) and loss rate (d). The image in plot (a) depicts the approximate posterior expectation for the entire time-evolution of the particle size density given by FIKS, whereas plots (e) and (f) illustrate the EKF and FIKS estimates corresponding to times 2 and 10 h, respectively. In plots (b)–(f), the blue and orange lines represent the approximate posterior expectations for EKF and FIKS, respectively, and the areas shaded with light blue and orange are the respective posterior error intervals. The true values of the corresponding quantities are drawn with a green line.

Figure 5 Cases 3 and 4 – steady-state system, showing true (simulated) processes: (a) particle size density (measurement), (b) nucleation rate, (c) growth rate and (d) wall loss rate.

In Case 1, the SNR is high, and – apart from the aforementioned exceptions – the estimate uncertainty is very low. Figure 4 shows the results of Case 2, where the SNR is significantly lower. As expected, the process rate estimates become less accurate than in Case 1. However, the change in the accuracy is quite small, especially for FIKS – demonstrating that the Kalman smoother estimates tolerate measurement noise rather well. As expected, the posterior error intervals of all estimated quantities are clearly wider than in the case of high SNR – this is a result of increased uncertainty in the particle counter observations. Both EKF and FIKS lead to safe posterior error intervals for the process rate parameters, and FIKS again gives clearly narrower posterior error intervals than EKF. The true values of the process rate parameters are once again within the posterior error intervals, further confirming the feasibility of Bayesian state estimation for quantifying the uncertainties of the process rates.

One may notice the appearance of spurious oscillations in the size distribution at low sizes during the particle formation event in Case 2 (Fig. 4a), which were completely smoothed out in Case 1 (Fig. 3a). These oscillations are due to the instability of Eqs. (5) and (6) which are corrected during the estimation by the assimilation of the data. However, while the data in Case 1 allow for the estimates to be neatly rectified, the data in Case 2 cannot completely make up for the instability. Note that the tolerance with respect to such modeling errors (e.g., discretization) can be further improved by so-called approximation error analysis (Huttunen and Kaipio, 2007).

3.2 Cases 3 and 4: steady state (SS)

3.2.1 SS: data simulation and parameter estimation

In the SS simulations in this paper, the condensation rate is both size dependent and time dependent, and the wall loss rate is time invariant but depends on size. Both the nucleation and condensation rates start from zero and grow within the first 30 min until they reach their stationary values. In such a case, the new particle formation and condensational growth are compensated for by wall losses, leading to the number concentration function to reach a steady state. The simulated process rates and the particle size density are illustrated in Fig. 5. Here, the size range of particles is [0.87 nm, 10.00 nm], and it is discretized into Q=1731 logarithmically distributed bins. Note that the value $d_{\mathrm{p}}^{*}=\mathrm{0.87}\mathrm{nm}$ is most certainly below any physically relevant critical sizes; however, from a simulation stand point, any size range could have been utilized to simulate a similar chamber experiment, only the values of the parameters would require adjustments.

Figure 6 Case 3 – SS, high SNR. State estimates for the particle size density (a, e, f), nucleation rate (b), growth rate (c) and wall loss rate (d). The images in plots (a) and (c) depict the approximate, FIKS-based posterior expectations for the entire time-evolution of the corresponding quantities. Plots (e) and (f) illustrate the EKF and FIKS estimates for the particle size density corresponding to times 1 and 2 h, respectively. In plots (b) and (d–f), the blue and orange lines represent the approximate posterior expectations for EKF and FIKS, respectively, and the areas shaded with light blue and orange are the respective posterior error intervals. The true values of the corresponding quantities are drawn with a green line.

Figure 7 Case 3 – SS and high SNR. Condensational growth rate estimates corresponding to four instants in time. The blue and orange lines represent the approximate posterior expectations for EKF and FIKS, respectively, and the areas shaded with light blue and orange are the respective posterior error intervals. The true growth rates are marked with green lines.

The synthetic SMPS data are simulated similarly to cases 1 and 2. Again, two cases corresponding to different SNRs are simulated by generating the Poisson-distributed observations corresponding to two sample volumes: V=200 cm³ (Case 3: SS, high SNR) and V=2 cm³ (Case 4: SS, low SNR). In Case 3, the range of SNR was [0,4440]; in Case 4 it was [0,44.4].

In the state space, the lower end of the size discretization bin (i.e., Ω₁ in Sect. 2.1.1) is centered at 1.1 nm (its lower boundary is about 1.08 nm). Many potential bins would separate the critical size $d_{\mathrm{p}}^{*}=\mathrm{0.87}\mathrm{nm}$ from the lower end of the state space $d_{\mathrm{p}}^{\text{min}}=\mathrm{1.08}\mathrm{nm}$. The reason that we are not setting $d_{\mathrm{p}}^{\text{min}}=d_{\mathrm{p}}^{*}$ is because, with the current technology, it is not realistic to assume that we can acquire reliable data down to the critical size. Note that it is possible, however, to extend the model down to the critical size in order to estimate the true nucleation rate – provided that the GDE is still valid for such small particles. However, without measurement from the lower end of the spectrum, the uncertainty will overwhelm the estimates, rendering them non-informative. This extension is outside of the scope of this paper.

The time evolution of the apparent particle formation and growth rates are modeled as second-order Markov processes (Eq. 12), whereas the wall loss rates are modeled as first-order Markov processes (Eq. 11). The growth and wall loss rates are assumed to be smooth with respect to size; hence, the state noises in the corresponding evolution models are modeled as correlated. For details of the model choices, we refer to Appendix B.

3.2.2 SS: results and discussion

Figure 6 illustrates the estimates of the particle size density and process rates in the case of high SNR (Case 3). For the particle size density and condensation rate, the respective images in Fig. 6a and c show the approximate FIKS-based posterior means of those time- and size-dependent variables. For the nucleation rate as well as for the instantaneous wall loss rate and size densities, both the filter and smoother estimates (approximate posterior expectations and error intervals) are plotted.

In this case, the smoother estimate for the particle size density is in very good correspondence with the true density (see Fig. 5a). Further, the uncertainty estimates for the particle size density are feasible: the true size density lies within the posterior error intervals given by EKF and FIKS.

The apparent particle formation rate is well estimated by both the EKF and the FIKS; the true value lies within the posterior error intervals most of the time. Only at the onset of the appearance of particles into the measured size range does the true value not lie within the uncertainty range, around t=25 min. Similarly to the NE cases, the FIKS estimates are less uncertain than those from the EKF, and the SNR levels are a major factor determining the accuracy of the estimates: the higher the SNR, the better the estimate. While the true nucleation rate plotted in dark green takes place at 0.87 nm, in the state-space model, the diameter that corresponds to the apparent particle formation, plotted in light green, is about 1.08 nm (the lower end of the smallest size bin, geometrical center at 1.1 nm). This difference between the model used for simulating the data and the model used in estimation causes the systematic error – underestimation of the nucleation rate.

Figure 8 Case 4 – SS and low SNR. State estimates for the particle size density (a, e, f), nucleation rate (b), growth rate (c) and wall loss rate (d). The images in plots (a) and (c) depict the approximate, FIKS-based posterior expectations for the entire time-evolution of the corresponding quantities. Plots (e) and (f) illustrate the EKF and FIKS estimates for the particle size density corresponding to times 1 and 2 h, respectively. In plots (b) and (d–f), the blue and orange lines represent the approximate posterior expectations for EKF and FIKS, respectively, and the areas shaded with light blue and orange are the respective posterior error intervals. The true values of the corresponding quantities are drawn with a green line.

Figure 9 Case 4 – SS and low SNR. Condensational growth rate estimates corresponding to four instants in time. The blue and orange lines represent the approximate posterior expectations for EKF and FIKS, respectively, and the areas shaded with light blue and orange are the respective posterior error intervals. The true growth rates are marked with green lines.

Figure 6c shows a clear trend in the quality of the smoother estimate for the condensational growth rate: at an early stage of the process (time ∼0.25 h), FIKS infers the growth rate reliably only in the smallest size classes (diameter ∼1 nm); in the larger particle sizes, the growth rate is heavily underestimated (see Fig. 5c). As time progresses, the growth rate estimates become gradually more reliable in larger and larger size classes.

The gradual improvement in the growth rate estimates in the larger size classes is a direct consequence of the propagation of the particle number density towards large size classes. Indeed, comparison of Fig. 6a and c reveals that the size class where FIKS catches the increase in the growth rate parameter (light/yellow area in the condensation rate image Fig. 6c) follows accurately the propagating front of the number density in Fig. 6a. The reason for this property of the growth rate estimate is obvious: in the size classes where the particle number density is very low, the measurement data do not carry information on the growth rate parameters. In the beginning of the process, the particle number density is low in all classes. When nucleation starts producing particles to the smallest size class and these particles grow, the growth is sensed by the particle size analyzer measurements. When the particle sizes keep increasing due to condensation, the measurements corresponding to increasingly larger size classes become sensitive to the process rate parameters.

Figure 7 illustrates the EKF and FIKS estimates of the growth rate corresponding to four instants in time. These plots confirm the above discussion on the growth rate estimation: at time 1 h and 6 min, both the EKF- and FIKS-based posterior means underestimate the growth rate in the large size classes (above ∼2.5 nm), and for the subsequent times, the EKF and especially FIKS estimates become reliable in gradually increasing size classes. Furthermore, Fig. 7 shows a trend in the evolution of the posterior error intervals of the growth rate: at time 1 h and 6 min, the posterior error intervals are really wide in classes >2.5 nm, reflecting high uncertainty in the growth rate estimates in this size range. As time progresses, the size range of low uncertainty spreads towards large size classes. This result demonstrates that Bayesian filtering and smoothing yield feasible posterior error estimates – indicating high uncertainty in the size ranges where the posterior expectations are unreliable in this example.

The EKF estimates of the growth rate in Fig. 7 show similar behavior to the FIKS estimates; the main differences are that, as expected, the posterior means given by EKF are more biased than those in FIKS, and the posterior error intervals are generally wider than with FIKS.

Figures 8 and 9 show the results of Case 4 – SS and low SNR data. The properties of the state estimates are very similar to those in Case 3, except for the anticipated differences: the lower accuracy of the approximate posterior means, and wider posterior error intervals – again indicating the increase in the uncertainty when the SNR of the measurements gets lower.

4 Conclusions

The radiative forcing caused by aerosols – and the underlying processes of new particle formation and growth – are currently understood to be the most uncertain factors in the prediction of climate change. While a significant effort has been directed into experimental campaigns for collecting data that carry indirect information on these processes and their effects, computational methods for analyzing the data are still insufficient for reliable estimation of the process rates.

In this paper, the problem of estimating the apparent particle formation, growth and loss rates of aerosols was cast in the framework of Bayesian state estimation. We modeled the dynamics of aerosol size distributions with the general dynamics equation and considered the (process rate) parameters in the model as unknown state variables. These size- and/or time-dependent variables were estimated together with the particle number density based on sequential particle counter measurements using the extended Kalman filter (EKF) and fixed interval Kalman smoother (FIKS). Furthermore, to quantify the uncertainties of the estimated variables, we also computed posterior error intervals for the process rate parameters.

The approach was tested with a set of numerical simulation studies, where two processes of different types were considered: (1) a nucleation event, which qualitatively represents a typical particle formation event in the atmosphere, and (2) a process approaching steady state, which is representative of a particle formation and growth in a chamber experiment. The EKF- and especially FIKS-based estimates for the process rates were generally reliable, even in cases of low SNR data. Moreover, the posterior error intervals were feasible – that is, the true process rates lie mostly within the approximate error intervals. The reason why FIKS estimates are generally superior to EKF estimates is that Bayesian smoothers also use future data when estimating a quantity at given time t, whereas filter estimates are based on only the data up to the point t. In online and control applications, of course, the future data are not available; however, in the applications targeted in this paper, the data analysis is done after the experiment – that is, the entire data set is available for Bayesian smoothing. The results of the numerical studies support the feasibility of the proposed approach to estimating the aerosol formation, growth and deposition rates, and quantifying their uncertainties. Based on these findings, we conclude that Bayesian state estimation (combined with aerosol particle dynamics modeling) offers a reliable tool for analyzing sequentially measured particle counter data. The estimated aerosol process rates and their uncertainties can improve the analysis of the experimental data and provide better insight on the particle formation and growth in the atmosphere. In the future, such analysis can potentially offer improved assessment of the radiative forcing by aerosols and its uncertainty. Eventually, this can lead to improved predictions and uncertainty quantification of the climate change via a combination of the enhanced process rate estimates and their uncertainties with global aerosol–climate models.

Appendix A: Parameters in data simulation

(Seinfeld and Pandis, 2016)(Seinfeld and Pandis, 2016)

Table A1 Summary of the models and parameters used for simulating data in nucleation event (NE) and steady-state (SS) cases. Here, the growth rate g is assumed to be a product of a size- and time-dependent parts – that is, g=g_d(d_p)g_t(d_t). “evol.” denotes “evolution”, and “obs.” represents “observed”.

Table A1 lists the models and parameters used for simulating the nucleation, condensation and deposition processes, and the particle size density, as well as modeling the particle counter data.

Appendix B: Model details and parameters in state estimation

Table B1 Qualitative description of the models used in state estimation for nucleation event (NE) and steady-state (SS) cases.

Figure B1 Covariance matrices of (a) the deposition rate λ in cases 1 and 2 and (b) growth rate g in cases 3 and 4.

In this Appendix, we describe the choices of models and parameters made for state estimation in the numerical examples (Sect. 3). We start by qualitatively summarizing the choices of models for the process rates in the two example cases, a nucleation event (NE) and steady state (SS), see Table B1.

Table B2 Parameters used in state estimation for nucleation event (NE) and steady-state (SS) cases.

In the following, we first describe how the smoothness of the size-dependent parameters is formulated. Next, we discuss modeling the time dependence of the parameters by first- and second-order Markov models. Finally, we list all parameter values chosen for each test case.

B1 Size-dependent processes: smoothness

The size-dependent process rate variables – loss rate (λ) in all test cases and the growth rate (g) in the SS cases – are assumed to be smooth functions of size. As these variables also depend on time, we model them as multivariate stochastic processes, particularly first-order Markov processes, as described below. The smoothness in size is accounted for by modeling elements of each of the associated random variables (size-discretized process rate variables ${\mathit{\xi }}_{\mathit{\phi }}^k=[{\mathit{\xi }}_{\mathit{\phi },\mathrm{1}}^k,\mathrm{\dots },{\mathit{\xi }}_{\mathit{\phi },Q}^k{]}^{\mathrm{T}}$, $\mathit{\phi }=\mathit{\lambda },g$) at a given time k as mutually correlated.

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{\mathbf{\Gamma }}_{{\mathit{\xi }}_{\mathit{\phi }}^k}(i,j)=& \sqrt{{\mathit{\sigma }}_{\mathit{\phi },i}^{\mathrm{2}}}\sqrt{{\mathit{\sigma }}_{\mathit{\phi },j}^{\mathrm{2}}}\exp \left(-{\displaystyle \frac{\mathrm{|}i-j\mathrm{|}}{{\mathit{\delta }}_{\mathit{\phi }}}}\right),\\ & \text{ for }\mathrm{1}\le i,j\le M\end{array}\end{array}$$

where ${\mathit{\sigma }}_{\mathit{\phi },i}^{\mathrm{2}}$ is the variance of ${\mathit{\xi }}_{\mathit{\phi },i}^k$, and δ_φ is a parameter defining how steeply the cross-covariance between elements ${\mathit{\xi }}_{\mathit{\phi },i}^k$ and ${\mathit{\xi }}_{\mathit{\phi },j}^k$ decreases as a function of the difference between indices i and j.

Figure B1 shows the covariance matrices of the two size-dependent variables: loss rate (λ) in cases 1 and 2 and growth rate (g) in cases 3 and 4. We note here that ${\mathit{\sigma }}_{\mathit{\lambda },i}^{\mathrm{2}}$ is constant in cases 1 and 2, whereas ${\mathit{\sigma }}_{g,i}^{\mathrm{2}}$ increases with particle size in cases 3 and 4. Details of parameter values are given in the following.

B2 Time-dependence: first-order Markov processes

Assume next that the process rate variable ${\mathit{\xi }}_{\mathit{\phi }}^k$ is modeled as a first-order Markov process. Here, we fix the covariance of ${\mathit{\xi }}_{\mathit{\phi }}^k$ first, as well as the temporal smoothness (first- or second-order Markov process). We then determine the covariance of the driving noise process.

When the covariance matrix ${\mathbf{\Gamma }}_{{\mathit{\xi }}_{\mathit{\phi }}^k}$ of ${\mathit{\xi }}_{\mathit{\phi }}^k$ is time invariant, we can calculate the covariance matrix of the state noise ${\mathit{\eta }}_{\mathit{\phi }}^k$ using Eq. (11) as follows:

$$\begin{array}{}\text{(B2)}& {\displaystyle }{\displaystyle }{\mathbf{\Gamma }}_{{\mathit{\eta }}_{\mathit{\phi }}^k}={\mathbf{\Gamma }}_{{\mathit{\xi }}_{\mathit{\phi }}^k}-{\mathbf{\Psi }}_{\mathit{\phi }}{\mathbf{\Gamma }}_{{\mathit{\xi }}_{\mathit{\phi }}^k}{\mathbf{\Psi }}_{\mathit{\phi }}^{\mathrm{T}}.\end{array}$$

In cases of size-dependent variables, we construct the covariance matrix ${\mathbf{\Gamma }}_{{\mathit{\xi }}_{\mathit{\phi }}^k}$ as described above, fix the state transition matrix Ψ_φ=r_φI, by choosing the parameter $r_{\mathit{\phi }}\in ]\mathrm{0},\mathrm{1}[$, and, finally, compute the state noise covariance matrix using Eq. (B2). Note, however, that in the case of a scalar variable modeled as a first-order Markov process φ, which does not depend on size, the state-transition matrix and the covariance are simply Ψ_φ=r_φ and ${\mathbf{\Gamma }}_{{\mathit{\xi }}_{\mathit{\phi }}^k}={\mathit{\sigma }}_{\mathit{\phi },i}^{\mathrm{2}}$, respectively.

B3 Time-dependence: second-order Markov processes

As shown in Table B1, the size-independent variables in cases 1 and 2 (g and J) are modeled as second-order Markov processes. For these variables, we choose the root parameters $r_{\mathit{\phi }}^{\mathrm{1}}$ and $r_{\mathit{\phi }}^{\mathrm{2}}$ in the evolution models using an approach adopted from the analysis of second-order systems, such as damped oscillators. We first define a characteristic time is T_φ and damping ratio ζ_φ, and then calculate $r_{\mathit{\phi }}^{\mathrm{1}}$ and $r_{\mathit{\phi }}^{\mathrm{2}}$ by solving the set of equations:

$$\begin{array}{}\text{(B3)}& {\displaystyle }& {\displaystyle }{r}_{\mathit{\phi }}^{\mathrm{1}}+r_{\mathit{\phi }}^{\mathrm{2}}=\mathrm{2}\left(\mathrm{1}-{\mathit{\zeta }}_{\mathit{\phi }}{\displaystyle \frac{\mathrm{2}\mathit{\pi }\mathrm{\Delta }{t}^k}{T_{\mathit{\phi }}}}\right)\text{(B4)}& {\displaystyle }& {\displaystyle }{r}_{\mathit{\phi }}^{\mathrm{1}}{r}_{\mathit{\phi }}^{\mathrm{2}}=\mathrm{1}-\mathrm{4}\mathit{\pi }{\mathit{\zeta }}_{\mathit{\phi }}{\displaystyle \frac{\mathrm{\Delta }{t}^k}{T_{\mathit{\phi }}}}+\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}{\left({\displaystyle \frac{\mathrm{\Delta }{t}^k}{T_{\mathit{\phi }}}}\right)}^{\mathrm{2}}.\end{array}$$

These parameters, along with separately chosen variances of the state noises ${\mathit{\eta }}_g^k$ and ${\mathit{\eta }}_J^k$ as well as the expectations and variances of the initial states ${\mathit{\xi }}_g^{\mathrm{0}}$ and ${\mathit{\xi }}_J^{\mathrm{0}}$, define the properties of the second-order Markov models.

B4 Parameter choices

All parameter values chosen for state estimation in each test case are listed in Table B2.