2.1.1. An idealized vertical prior

In passive satellite sensors, AOD retrieval algorithms are used to produce an estimate of the AOD from measurements of top-of-atmosphere reflectance. They are typically based on a forward model of top-of-atmosphere reflectance (Wu et al., Reference Wu, de Graaf and Menenti2017). The parameters required to run this forward model are obtained using a look-up-table algorithm, which admits as indices aerosol properties such as the aerosol type, the AOD, and the form of the vertical profile. When the satellite sensor captures a measure of reflectance, the look-up-table can be used to solve the corresponding inverse problem and estimate the AOD, provided the aerosol properties have been predetermined. Consequently, to address the inverse problem, AOD retrieval algorithms need to assume a form for aerosol vertical profile.

The form of the vertical profile is typically idealized, using either a Gaussian profile, or in the simplest case an exponential profile (Wu et al., Reference Wu, de Graaf and Menenti2017; Li et al., Reference Li, Li, Dubovik, Zeng and Yung2020). The exponential profile takes the form $ {b}_{\mathrm{ext}}(h)\hskip0.22em \propto \hskip0.22em {e}^{-h/L} $ , where $ L $ is a fixed height scale parameter that is typically taken as the top altitude of the boundary layer ( $ 2\hskip0.22em \mathrm{km} $ )—although different algorithms may adopt varying values.

Whilst idealized, these profiles capture a key element of aerosol vertical distribution: most CCN lies at low altitude in the boundary layer. This constitutes a form of domain knowledge, which we argue is important to encode within the design of a prior over aerosol vertical profiles. Thus, we propose to use an idealized exponential profile as a structural component of the prior placed over the extinction profile.

2.1.2. Weighting the ideal profile

Clouds’ local meteorology influences clouds properties, and hence the sign and magnitude of the effective radiative forcing due to ACIs (Ackerman et al., Reference Ackerman, Kirkpatrick, Stevens and Toon2004; Small et al., Reference Small, Chuang, Feingold and Jiang2009; Chen et al., Reference Chen, Christensen, Stephens and Seinfeld2014; Douglas and L’Ecuyer, Reference Douglas and L’Ecuyer2020)—this source of heterogeneity is at the heart of the uncertainty over how clouds impact climate.

The impact of local meteorology on ACIs can be characterized by the following set of environmental confounders (Köhler, Reference Köhler1936; Twomey, Reference Twomey1974): temperature $ T $ , pressure $ P $ , relative humidity RH, and vertical velocity or updraft $ \omega $ . Standing as a proxy for aerosols, the AOD is also impacted by its surrounding meteorology (Christensen et al., Reference Christensen, Neubauer, Poulsen, Thomas, McGarragh, Povey, Proud and Grainger2017; Jesson et al., Reference Jesson, Douglas, Manshausen, Meinshausen, Stier, Gal and Shalit2022). These meteorological variables should thus also modulate the extinction coefficient.

A great advantage of these simple meteorological variables is that they (or their proxies) are readily available on multiple pressure levels in reanalysis data. Reanalysis data, such as ERA5 (Muñoz-Sabater et al., Reference Muñoz-Sabater, Dutra, Agustí-Panareda, Albergel, Arduini, Balsamo, Boussetta, Choulga, Harrigan, Hersbach, Martens, Miralles, Piles, Rodríguez-Fernández, Zsoter, Buontempo and Thépaut2021), combines observations and model simulations through physical laws to provide the most accurate representation of past and present climate and meteorology. Hence, $ T $ , $ P $ , RH, and $ \omega $ can be reliably used as vertically resolved predictors. Finally, note that $ {b}_{\mathrm{ext}} $ is also a spatiotemporal field which exhibits smoothness across spatial and temporal dimensions. Such regularity should be included in the modeling.

Let $ x\mid h $ denote a $ d $ -dimensional vector resulting from the concatenation of spatiotemporal and meteorological variables for a given altitude $ h $ . For example, one can consider

(2) $$ x=\left(t,\mathrm{lat},\mathrm{lon},T,P,\mathrm{RH},\omega \right), $$

where lat and lon respectively denote latitude and longitude. In the remainder of this paper, we will denote $ \mathcal{X}\hskip0.22em \subseteq \hskip0.22em {\mathrm{\mathbb{R}}}^d $ the space in which $ x $ takes values.

We propose to model the extinction coefficient $ {b}_{\mathrm{ext}} $ by weighting an idealized exponential vertical profile with a positive weight function $ w:\mathcal{X}\to \left(0,+\infty \right) $ . Namely, the chosen prior for the extinction coefficient profile is denoted $ \varphi $ and takes the simple form

(3) $$ \varphi \left(x|h\right)=w\left(x|h\right){e}^{-h/L}. $$

This weight function is meant to capture finer details of variability in the extinction profile, putting more mass in regions where meteorological predictors suggest aerosol loading is likely to be higher.

2.1.3. Probabilistic modeling of the weighting function

Whilst the extinction coefficient should indeed be modulated by $ x\mid h $ , we expect this relationship to be non-trivial and highly non-linear (Jesson et al., Reference Jesson, Douglas, Manshausen, Meinshausen, Stier, Gal and Shalit2022). For this reason, we propose to learn the weighting function $ w $ using non-linear (e.g., kernel-based) statistical machine-learning methodologies. Furthermore, we want the prior to reflect our lack of knowledge about the relationship between $ x\mid h $ and $ {b}_{\mathrm{ext}}(h) $ . To account for this epistemic uncertainty, we propose a Bayesian formulation of the weighting function.

GPs (Rasmussen and Williams, Reference Rasmussen and Williams2005) are a ubiquitous class of expressive Bayesian priors over real-valued functions. They have been widely used in various nonlinear and nonparametric regression problems in geosciences (Camps-Valls et al., Reference Camps-Valls, Verrelst, Muñoz-Marí, Laparra, Mateo-Jimenez and Gómez-Dans2016). A GP is fully determined by its mean function and its covariance function. The covariance function—called kernel—is typically user-specified as a positive definite bivariate function on the input data.

We place a GP prior over the weight function. To ensure the weights are strictly positive, we further warp the GP with a positive transform $ \psi :\mathrm{\mathbb{R}}\to \left(0,+\infty \right) $ . Formally, let $ m:\mathcal{X}\to \mathrm{\mathbb{R}} $ be a mean function and $ k:\mathcal{X}\times \mathcal{X}\to \mathrm{\mathbb{R}} $ denote a positive definite kernel, we model the weighting function as

(4) $$ w\left(x|h\right)=\psi \left(f\left(x|h\right)\right)\hskip1em \mathrm{where}\hskip1em f\sim \mathrm{GP}\left(m,k\right). $$

The Bayesian prior defined by $ f $ represents a latent variable that maps through the positive link function $ \psi $ onto the weight function. A simple choice for $ \psi $ is the exponential function, which as we will see in Section 3 is a mathematically convenient choice.

Whilst $ \psi \hskip0.35em \circ f $ describes an expressive probability distribution over nonlinear positive functions, it remains interpretable. The choice of kernel specifies how covariance is encoded, that is, $ \mathrm{Cov}\left(f(x),f\left({x}^{\prime}\right)\right)=k\left(x,{x}^{\prime}\right) $ and controls the functional class the GP belongs to. For example, the Matérn class of covariance functions offers control over the functional smoothness of the GP (Stein, Reference Stein1999; Rasmussen and Williams, Reference Rasmussen and Williams2005).

2.2.1. Relevance of the log-normal distribution

The log-normal distribution has been reported to provide a good fit to the AOD in studies focusing on locations in North America and Europe (O’Neill et al., Reference O’Neill, Ignatov, Holben and Eck2000). This distribution is particularly suitable because the AOD is strictly positive and exhibits a right-skewed distribution concentrated around small values ( $ \approx $ 0.14).

The log-normal distribution is a right-skewed continuous probability distribution with support $ \left(0,+\infty \right) $ . It is specified by two parameters: a location parameter $ \mu \in \mathrm{\mathbb{R}} $ , and a scale parameter $ \sigma >0 $ . In particular, if $ Z\sim \mathcal{N}\left(0,1\right) $ is a standard normal random variable, then $ {e}^{\mu +\sigma Z} $ follows a log-normal distribution with location $ \mu $ and scale $ \sigma $ , which we denote $ \mathrm{\mathcal{L}}\mathcal{N}\left(\mu, \sigma \right) $ .

It is often preferable to express the log-normal distribution as a function of its mean rather than its location, as it is more interpretable and offers a convenient parametrization for including higher-order model parameters. Formally, let $ \eta =\unicode{x1D53C}\left[\mathrm{\mathcal{L}}\mathcal{N}\left(\mu, \sigma \right)\right] $ denote the expectation of the distribution and suppose the scale $ \sigma $ is known. The $ \eta $ can be expressed using $ \mu $ and $ \sigma $ as $ \eta ={e}^{\mu +{\sigma}^2/2}\hskip0.35em \iff \hskip0.35em \mu =\log \eta -\frac{\sigma^2}{2} $ , which yields the mean-reparametrized formulation of the log-normal distribution

(5) $$ \mathrm{\mathcal{L}}\mathcal{N}\left(\mu, \sigma \right)=\mathrm{\mathcal{L}}\mathcal{N}\left(\log \eta -\frac{\sigma^2}{2},\sigma \right). $$

We collect AERONET sun-photometers AOD measurements from 1315 stations between 1993 and 2021 (Holben et al., Reference Holben, Eck, Slutsker, Tanré, Buis, Setzer, Vermote, Reagan, Kaufman, Nakajima, Lavenu, Jankowiak and Smirnov1998), and shown in Figure 1 that the empirical marginal distribution of AERONET observations can indeed be appropriately fitted using a log-normal density. Whilst Figure 1 only depicts a marginal distribution for the AOD, we draw motivation from this to formulate a log-normal observation model for $ {b}_{\mathrm{ext}} $ and $ \tau $ conditional on $ \varphi \left(x|h\right) $ .

Figure 1. Left: Log-normal density (red) fitted with maximum likelihood estimates to the empirical distribution of AERONET AOD at 500 nm ( $ {\tau}_{500} $ ) from 1315 stations between 1993 and 2021 (green). Right: logspace plot of the left panel. It demonstrates a sound fit for the normal distribution in the logspace, albeit with a slight right skew; $ \hat{\mu}=-2.06 $ , $ \hat{\sigma}=0.96 $ .

2.2.2. An observation model for b_ext

Whilst we do not usually observe $ {b}_{\mathrm{ext}}(h) $ it is still useful to specify an observation model for the extinction coefficient. Indeed, at inference stage, $ \varphi \left(x|h\right) $ might fail to capture observational noise which will inevitably hinder variance calibration for the predicted extinction coefficient profile distribution. Using an observation model over $ {b}_{\mathrm{ext}} $ introduces an additional degree of freedom that will help calibrate the variance of observations.

We propose a simple log-normal observation model with mean $ \varphi \left(x|h\right) $ for the extinction coefficient,

(6) $$ {b}_{\mathrm{ext}}(h)\mid \varphi \left(x|h\right)\sim \mathrm{\mathcal{L}}\mathcal{N}\left(\log \varphi \left(x|h\right)-\frac{\sigma_{\mathrm{ext}}^2}{2},{\sigma}_{\mathrm{ext}}\right), $$

where $ {\sigma}_{\mathrm{ext}}>0 $ is a scale parameter. This observation model conserves the meanFootnote ² from our prior $ \varphi \left(x|h\right) $ and simply includes observational noise through the scale parameter $ {\sigma}_{\mathrm{ext}} $ . When extinction coefficient observations are available, $ {\sigma}_{\mathrm{ext}} $ can be calibrated following the procedure described in Section 3.3.3.

We emphasize that tuning $ {\sigma}_{\mathrm{ext}} $ requires only a small number of extinction coefficient observations. This is in contrast to the AOD observations, which constitute the principal source of data we observe in the context of this work.

2.2.3. An observation model for τ

AOD observations constitute the primary source of information on aerosols optical properties in the context of this work since it is collected routinely at a global scale by satellite products. Satellite retrievals are however prone to distortion, that may arise from observational noise (Remer et al., Reference Remer, Kaufman, Tanré, Mattoo, Chu, Martins, Li, Ichoku, Levy, Kleidman, Eck, Vermote and Holben2005) or assumptions made in the AOD retrieval algorithm (Mielonen et al., Reference Mielonen, Levy, Aaltonen, Komppula, de Leeuw, Huttunen, Lihavainen, Kolmonen, Lehtinen and Arola2011; Wu et al., Reference Wu, De Graaf and Menenti2016). Therefore, as for the extinction coefficient, we propose to model these distortions using a log-normal observation model for the AOD conditional on $ \varphi \left(x|h\right) $ .

Formally, let $ \eta ={\int}_0^H\varphi \left(x|h\right)\hskip0.1em \mathrm{d}h $ be the column integration of our prior, then we propose a mean-parametrized observation model for the AOD conditional on $ \eta $ , given by

(7) $$ \tau \mid \eta \sim \mathrm{\mathcal{L}}\mathcal{N}\left(\log \eta -\frac{\sigma^2}{2},\sigma \right), $$

(8) $$ \eta ={\int}_0^H\varphi \left(x|h\right)\hskip0.1em \mathrm{d}h, $$

with a scale parameter $ \sigma >0 $ which can be tuned against AOD observations jointly with other model parameters following the procedure described in Section 2.3.

When working with multiple AOD observations $ {\tau}_1,\dots, {\tau}_n $ , the scale parameter $ \sigma $ will be assumed shared among atmospheric columns. The mean parameter $ \eta $ (or location parameter $ \mu $ for the canonical parametrization) will however be column-specific.

2.3.1. Finite-sample problem formulation

We now make things concrete and assume we observe the AOD for $ n $ columns, which we stack into the vector $ \tau ={\left[{\tau}_1\hskip0.5em \dots \hskip0.5em {\tau}_n\right]}^{\top}\in {\mathrm{\mathbb{R}}}^n $ . For the $ i $ th column, we also observe $ {m}_i $ vertically resolved meteorological covariates $ {x}_i^{(1)},\dots, {x}_i^{\left({m}_i\right)} $ and their respective altitudes $ {h}_i^{(1)}<\dots <{h}_i^{\left({m}_i\right)} $ , such that $ {x}_i^{(j)}\sim p\left(x|{h}_i^{(j)}\right) $ . We concatenate these observations into a dataset $ \mathcal{D}={\left\{{\left({x}_i^{(j)},{h}_i^{(j)}\right)}_{j=1}^{m_i},{\tau}_i\right\}}_{i=1}^n $ and denote $ M={\sum}_{i=1}^n{m}_i $ the total number of vertically resolved samples. The model description for the $ i $ th column is summarized in Figure 2.

Figure 2. Observation models and prior formulation for the $ i $ ^th atmospheric column. The model follows a hierarchical Bayesian structure represented by the arrows. We first specify the prior, and then express the observation models conditionally on the prior.

Our objective is to use $ \mathcal{D} $ to learn the mapping $ \varphi \left(x|h\right) $ . Within the Bayesian framework, this objective is twofold:

1. (A) Update the prior placed on $ \varphi $ with the observations from $ \mathcal{D} $ . This corresponds to computing the posterior distribution of $ \varphi $ given $ \tau $ .

2. (B) Tune the model hyperparameters by maximizing the marginal log-likelihood $ \log p\left(\tau \right) $ . The hyperparameters include the log-normal scale parameter $ \sigma $ and any parameter from the GP mean m and kernel $ k $ .

Regarding point (A), since $ \varphi $ results from a transformation of GP $ f $ , we will rather focus on the posterior distribution of the GP directly for convenience. If $ \mathrm{f}\in \mathrm{\mathbb{R}} $ denotes a realization of $ f $ at any input $ x\mid h $ , the posterior distribution of $ f\left(x|h\right) $ given observations is denoted $ p\left(\mathrm{f}|\tau \right) $ . Having access to the posterior $ p\left(\mathrm{f}|\tau \right) $ allows computation of the predictive mean and variance of $ \varphi \left(x|h\right) $ following

(9) $$ \unicode{x1D53C}\left[\varphi \left(x|h\right)|\tau \right]={\int}_{\mathrm{\mathbb{R}}}\psi \left(\mathrm{f}\right){e}^{-h/L}p\left(\mathrm{f}|\tau \right)\hskip0.1em \mathrm{df}, $$

(10) $$ \mathrm{Var}\left(\varphi \left(x|h\right)|\tau \right)=\unicode{x1D53C}\left[\varphi {\left(x|h\right)}^2|\tau \right]-\unicode{x1D53C}{\left[\varphi \left(x|h\right)|\tau \right]}^2. $$

The above can be estimated with Monte Carlo by drawing samples from the posterior $ p\left(\mathrm{f}|\tau \right) $ . In Section 3.3.1, we will obtain a simpler closed-form solution for the particular choice $ \psi =\exp $ .

Regarding point (B), the marginal log-likelihood $ \log p\left(\tau \right) $ is unfortunately intractable. Indeed, let $ \mathbf{x}={\left[\begin{array}{ccc}{x}_1^{(1)}& \dots & {x}_n^{\left({m}_n\right)}\end{array}\right]}^{\top}\in {\mathcal{X}}^M $ denote the concatenation of all input entries from the dataset and let $ \mathbf{f}\in {\mathrm{\mathbb{R}}}^M $ denote a realization of $ f\left(\mathbf{x}\right) $ . Assuming independence of AOD observations $ {\tau}_1,\dots, {\tau}_n $ conditionally on the GP realization over atmospheric columns $ \mathbf{f} $ , the marginal likelihood $ p\left(\tau \right) $ is expressed in terms of the observation model and prior distributions following

(11) $$ p\left(\tau \right)={\int}_{{\mathrm{\mathbb{R}}}^M}p\left(\tau |\mathbf{f}\right)p\left(\mathbf{f}\right)\hskip0.1em \mathrm{d}\mathbf{f}\hskip1em \mathrm{with}\hskip1em p\left(\tau |\mathbf{f}\right)=\prod \limits_{i=1}^np\left({\tau}_i|{\mathbf{f}}_i\right), $$

where $ {\mathbf{f}}_i={\left[\begin{array}{ccc}{\mathbf{f}}_i^{(1)}& \dots & {\mathbf{f}}_i^{\left({m}_i\right)}\end{array}\right]}^{\top}\in {\mathrm{\mathbb{R}}}^{m_i} $ corresponds to the GP realization over the $ i $ th column only. Because the observation model is log-normal, this integral is however intractable. To circumvent this intractability, we propose to tackle objective (B) by maximizing a proxy of the marginal log-likelihood presented in Section 2.3.3.

2.3.2. A sparse variational approximation of the posterior

We start by addressing objective (A). The predictive posterior distribution of interest is given by

(12) $$ p\left(\mathrm{f}|\tau \right)=\frac{p\left(\tau |\mathrm{f}\right)p\left(\mathrm{f}\right)}{\int_{\mathrm{\mathbb{R}}}p\left(\tau |\mathrm{f}\right)p\left(\mathrm{f}\right)\hskip0.1em \mathrm{df}}. $$

The integral denominator is not available in closed form, making the posterior intractable. We propose to use a variational approximation scheme (Titsias, Reference Titsias2009; Matthews et al., Reference Matthews, Hensman, Turner and Ghahramani2016; Leibfried et al., Reference Leibfried, Dutordoir, John and Durrande2020) to substitute this intractable inference problem with a tractable optimization problem. In addition, we make the variational approximation sparse (Titsias, Reference Titsias2009) such that the model can scale to large amounts of data.

Let $ \mathbf{z}={\left[\begin{array}{ccc}{z}_1& \dots & {z}_p\end{array}\right]}^{\top}\in {\mathcal{X}}^p $ be a set of $ p\ll M $ inducing locations over the space of inputs. Their evaluation by the GP follows a multivariate normal distribution $ f\left(\mathbf{z}\right)\sim \mathcal{N}\left(0,{\mathbf{K}}_{\mathbf{zz}}\right) $ , where $ {\mathbf{K}}_{\mathbf{zz}}=k\left(\mathbf{z},\mathbf{z}\right) $ . We denote $ \mathbf{u}=f\left(\mathbf{z}\right)\in {\mathrm{\mathbb{R}}}^p $ and refer to this vector as inducing variables.

A $ p $ -dimensional parametric distribution is set over these inducing variables. We choose this distribution as a multivariate normal defined by $ q\left(\mathbf{u}\right)\hskip0.35em := \hskip0.35em \mathcal{N}\left(\mathbf{u}|{\boldsymbol{\mu}}_{\mathbf{z}},{\boldsymbol{\Sigma}}_{\mathbf{z}}\right) $ . $ {\boldsymbol{\mu}}_{\mathbf{z}}\in {\mathrm{\mathbb{R}}}^p,{\boldsymbol{\Sigma}}_{\mathbf{z}}\in {\mathrm{\mathbb{R}}}^{p\times p} $ are called the variational parameters and need to be tuned such that $ q\left(\mathbf{u}\right) $ best approximates the true posterior $ p\left(\mathbf{u}|\tau \right) $ .

Once this is achieved, we take as an approximation to $ p\left(\mathrm{f}|\tau \right) $ the variational posterior defined by $ q\left(\mathrm{f}\right)\hskip0.35em := \hskip0.35em {\int}_{{\mathrm{\mathbb{R}}}^p}p\left(\mathrm{f}|\mathbf{u}\right)p\left(\mathbf{u}\right)\hskip0.1em \mathrm{d}\mathbf{u} $ , which is given in closed form by

(13) $$ q\left(\mathrm{f}\right)=\mathcal{N}\left(\mathrm{f}|{\overline{\mu}}_{x\mid h},{\overline{\Sigma}}_{x\mid h}\right), $$

(14) $$ {\overline{\mu}}_{x\mid h}=k\left(x|h,\mathbf{w}\right){\mathbf{K}}_{\mathbf{z}\mathbf{z}}^{-1}{\boldsymbol{\mu}}_{\mathbf{z}}, $$

(15) $$ {\overline{\Sigma}}_{x\mid h}=k\left(x|h,x|h\right)-k\left(x|h,\mathbf{z}\right)\left({\mathbf{K}}_{\mathbf{z}\mathbf{z}}^{-1}-{\mathbf{K}}_{\mathbf{z}\mathbf{z}}^{-1}{\boldsymbol{\Sigma}}_{\mathbf{z}}{\mathbf{K}}_{\mathbf{z}\mathbf{z}}^{-1}\right)k\left(\mathbf{z},x|h\right). $$

Naturally, (13) can be extended to describe a variational posterior over multiple GP entries. Namely, let $ {\mathbf{x}}^{\ast}\in {\mathcal{X}}^D $ denote a vector of input entries. If $ {\mathbf{f}}^{\ast}\in {\mathrm{\mathbb{R}}}^D $ denotes a realization of $ f\left({\mathbf{x}}^{\ast}\right) $ , then the associated variational posterior is given by

(16) $$ q\left({\mathbf{f}}^{\ast}\right)=\mathcal{N}\left({\mathbf{f}}^{\ast }|{\overline{\mu}}_{{\mathbf{x}}^{\ast }},{\overline{\boldsymbol{\Sigma}}}_{{\mathbf{x}}^{\ast }}\right), $$

(17) $$ {\overline{\mu}}_{{\mathbf{x}}^{\ast }}=k\left({\mathbf{x}}^{\ast },\mathbf{w}\right){\mathbf{K}}_{\mathbf{z}\mathbf{z}}^{-1}{\mu}_{\mathbf{z}}\in {\mathrm{\mathbb{R}}}^D, $$

(18) $$ {\overline{\boldsymbol{\Sigma}}}_{{\mathbf{x}}^{\ast }}=k\left({\mathbf{x}}^{\ast },{\mathbf{x}}^{\ast}\right)-k\left({\mathbf{x}}^{\ast },\mathbf{z}\right)\left({\mathbf{K}}_{\mathbf{z}\mathbf{z}}^{-1}-{\mathbf{K}}_{\mathbf{z}\mathbf{z}}^{-1}{\boldsymbol{\Sigma}}_{\mathbf{z}}{\mathbf{K}}_{\mathbf{z}\mathbf{z}}^{-1}\right)k\left(\mathbf{z},{\mathbf{x}}^{\ast}\right)\in {\mathrm{\mathbb{R}}}^{D\times D}. $$

The sparse nature of this approach becomes apparent in (17) and (18). Indeed, regardless of the number of samples we wish to evaluate the variational posterior over, we only need to invert a $ p\times p $ matrix, incurring a $ \mathcal{O}\left({p}^3\right) $ computational cost.

2.3.3. Learning the variational parameters $ {\mu}_{\mathbf{z}},{\boldsymbol{\Sigma}}_{\mathbf{z}} $

As mentioned above, the variational parameters $ {\boldsymbol{\mu}}_{\mathbf{z}},{\boldsymbol{\Sigma}}_{\mathbf{z}} $ need to be tuned such that $ q\left(\mathbf{u}\right) $ best approximates the posterior $ p\left(\mathbf{u}|\tau \right) $ , which is intractable. This problem is casted as the maximization of an objective called the evidence lower-bound (ELBO), given by

(19) $$ \mathrm{ELBO}(q)={\unicode{x1D53C}}_{q\left(\mathbf{f}\right)}\left[\log p\left(\tau |\mathbf{f}\right)\right]-\mathrm{KL}\left(q\left(\mathbf{u}\right)\Big\Vert p\left(\mathbf{u}\right)\right). $$

The ELBO is a lower-bound to the marginal log-likelihood $ \log p\left(\tau \right)\ge \mathrm{ELBO}(q) $ . It can thus be used as a proxy of $ \log p\left(\tau \right) $ to also tune the model hyperparameters and fulfill objective (B).

The second term in (19) is the Kullback–Leibler divergence between two multivariate normal distributions. It admits a closed-form expression and can be computed. The first term, on the other hand, is an expected log-likelihood under the variational posterior which cannot be analytically computed. It can be decomposed into column-wise terms

(20) $$ {\unicode{x1D53C}}_{q\left(\mathbf{f}\right)}\left[\log p\left(\tau |\mathbf{f}\right)\right]=\sum \limits_{i=1}^n{\unicode{x1D53C}}_{q\left({\mathbf{f}}_i\right)}\left[\log \mathrm{\mathcal{L}}\mathcal{N}\left({\tau}_i\left|\hskip0.1em \log {\eta}_i-\frac{\sigma^2}{2},\sigma \right.\right)\right]. $$

To estimate (20), we must first evaluate the mean parameter $ {\eta}_i={\int}_0^H\psi \left(f\left({x}_i|h\right)\right){e}^{-h/L}\hskip0.1em \mathrm{d}h $ with the finite number of GP evaluations $ {\mathbf{f}}_i $ we have access to. We propose to use the trapezoidal integration scheme given by

(21) $$ {\hat{\eta}}_i=\sum \limits_{j=1}^{m_i-1}\frac{\psi \left({\mathbf{f}}_i^{\left(j+1\right)}\right){e}^{-{h}_i^{\left(j+1\right)}/L}-\psi \left({\mathbf{f}}_i^{(j)}\right){e}^{-{h}_i^{(j)}/L}}{2}\left({h}_i^{\left(j+1\right)}-{h}_i^{(j)}\right). $$

While we choose the trapezoidal rule for simplicity, we note that alternative finite-sample integration schemes can be chosen here in accordance with the needs.

Second, because of the log-normal observation model, the expected log-likelihood remains intractable. To estimate it, while allowing backpropagation through the variational parameters, we use a reparametrization trick (Kingma et al., Reference Kingma, Salimans and Welling2015). Namely, we sample $ {\boldsymbol{\unicode{x025B}}}_i\sim \mathcal{N}\left(0,{\mathbf{I}}_{m_i}\right) $ and compute $ {\mathbf{f}}_i={\overline{\boldsymbol{\mu}}}_i+{\overline{\boldsymbol{\Sigma}}}_i^{1/2}{\boldsymbol{\unicode{x025B}}}_i $ , where $ {\overline{\boldsymbol{\mu}}}_i,{\overline{\boldsymbol{\Sigma}}}_i $ are the variational posterior parameters for the $ i $ th column and are obtained by application of (17) and (18) over the predictors of the $ i $ th column. The resulting GP sample $ {\mathbf{f}}_i $ is then used to estimate the mean parameter $ {\hat{\eta}}_i $ following (21) and we can approximate the expected log-likelihood with its one-sample estimate

(22) $$ {\unicode{x1D53C}}_{q\left({\mathbf{f}}_i\right)}\left[\log \mathrm{\mathcal{L}}\mathcal{N}\left({\tau}_i\left|\hskip0.1em \log {\eta}_i-\frac{\sigma^2}{2},\sigma \right.\right)\right]\approx \log \mathrm{\mathcal{L}}\mathcal{N}\left({\tau}_i\left|\hskip0.1em \log {\hat{\eta}}_i-\frac{\sigma^2}{2},\sigma \right.\right). $$

This method allows estimation of the ELBO objective, which in turn can be maximized with respect to the variational parameters $ {\boldsymbol{\mu}}_{\mathbf{z}},{\boldsymbol{\Sigma}}_{\mathbf{z}} $ using a stochastic gradient approach for example. As mentioned above, the model hyperparameters can also be tuned jointly with this objective, hence fulfilling objective (B). These include the log-normal scale $ \sigma $ or the kernel $ k $ hyperparameters (e.g., variances and lengthscales), with an option to parametrize these kernels using feature maps given by deep neural networks (Law et al., Reference Law, Zhao, Chan, Huang and Sejdinovic2019). As it is standard in sparse variational GPs (Titsias, Reference Titsias2009), we also learn the inducing locations $ \mathbf{z} $ .

3.1.1. The ECHAM-HAM dataset

The aerosol-climate model ECHAM-HAM is a self-consistent global climate model of aerosol radiative properties and CCN which demonstrates excellent agreement with AOD measurements from ground-based sun-photometers and satellite retrievals (Stier et al., Reference Stier, Feichter, Kinne, Kloster, Vignati, Wilson, Ganzeveld, Tegen, Werner, Balkanski, Schulz, Boucher, Minikin and Petzold2005). It computes the evolution of log-normal aerosol mass and number modes—for species sulfates, black carbon, organic carbon, sea salt, and mineral dust—by taking into account physical and chemical particle processes. The simulation used includes aerosol optical properties, aerosol tracers, and meteorological variables at a 1.8°  $ \times $  1.8° horizontal resolution on a Gaussian grid and over 31 levels of vertical resolution and for 8 regularly spaced time steps over a day (06/06/2008). The resulting dataset counts $ 147,456 $ atmospheric columns used as training points. Table 1 provides its detailed dimensions.

Table 1. Gridded variables from ECHAM-HAM simulation data

Note. The grid includes 8 time steps ( $ t $ ), 96 latitude levels (lat), 192 longitude levels (lon), and 31 vertical pressure levels (lev)—which is a proxy for $ h $ . Our objective is to vertically disaggregate the response $ \tau $ using the vertically resolved predictors ( $ T,P,\mathrm{RH},\omega $ ) and spatiotemporal columns locations ( $ t $ , lat, lon).

The dataset captures two important aspects of aerosol properties variability: (i) because the dataset is global, it captures the spatial variability of aerosols, which is arguably the greatest contribution to the overall variability of aerosol properties. Indeed, the distribution of particles in the troposphere is highly heterogeneous due to the combination of spatially diverse sources and the relatively short lifetime of particles and gases (Carslaw and Pringle, Reference Carslaw, Pringle and Carslaw2022); (ii) the dataset covers a full day, which provides a good coverage of the diurnal cycle of aerosols.

However, a single day of data may fail to account for the seasonal variations of aerosol properties. The seasonal variability of aerosols is strongly regionally dependent, where variations in aerosol concentration can exceed a factor of 10 in regions prone to episodic emissions (e.g., wildfires, human activity), but vary only by 10% over a year in regions with persistent atmospheric circulation such as the Pacific subtropics (Carslaw and Pringle, Reference Carslaw, Pringle and Carslaw2022). Therefore, we expect that the demonstrated predictive performance of our model over this dataset should carry over to more stable regions for different seasons, but it remains uncertain from experiments on this dataset whether the model may generalize to the regional temporal variability induced by episodic events.

3.1.2. Choice of variables

ECHAM-HAM offers simulations of aerosol properties at particular wavelengths that match the wavelengths at which operate specific instruments for evaluation. We use the ECHAM-HAM AOD at 550 nm as the response variable to vertically disaggregate and the ECHAM-HAM extinction coefficient at 533 nm as the groundtruth variable to evaluate our predictions against. The AOD and the extinction coefficients are simulated in ECHAM-HAM at slightly different wavelengths in order to match the wavelength at which operate the respective sensors (550 nm for the MODIS AOD product; Remer et al., Reference Remer, Kaufman, Tanré, Mattoo, Chu, Martins, Li, Ichoku, Levy, Kleidman, Eck, Vermote and Holben2005 and 532 nm for the CALIOP extinction; Winker et al., Reference Winker, Tackett, Getzewich, Liu, Vaughan and Rogers2013).

In order to make the model easily applicable to any satellite AOD retrieval, we select vertically resolved predictors among standard meteorological variables, that could in practice be easily obtained from reanalysis data: $ T $ , $ P $ , RH, and $ \omega $ . By limiting ourselves to these variables, which are comprehensively covered by reanalysis products, we favor a global applicability of the model to AOD observations from most times and locations.

In what follows, the input variable writes $ x=\left(t,\mathrm{lat},\mathrm{lon},T,P,\mathrm{RH},\omega \right) $ . We standardize meteorological predictors and transform them using a rank-based inverse normal mapping to achieve approximately normal distributions.

3.2.1. Choice of idealized profile height scale L

A simple choice for the idealized profile height scale is the boundary layer upper altitude, typically taken at $ 2\hskip0.22em \mathrm{km} $ . This is the value assumed by MISR, VIIRS, and MODIS C6_DT AOD retrieval algorithms (Levy et al., Reference Levy, Remer, Mattoo, Vermote and Kaufman2007; Kahn and Gaitley, Reference Kahn and Gaitley2015; Laszlo and Liu, Reference Laszlo and Liu2016). In our experiment, we also choose to set $ L=2\hskip0.22em \mathrm{km} $ .

3.2.2. Choice of mean and covariance structure

We set the GP with a zero mean function $ m=0 $ and kernel taken as the sum of a spatiotemporal and a meteorological contribution to the covariance

(23) $$ k\left(x|h,{x}^{\prime }|{h}^{\prime}\right)={\gamma}_1{k}_{\mathrm{ST}}\left(\left\{t,\mathrm{la}\mathrm{t},\mathrm{lo}\mathrm{n}\right\},\Big\{{t}^{\prime },\mathrm{la}{\mathrm{t}}^{\prime },\mathrm{lo}{\mathrm{n}}^{\prime}\Big\}\right)+{\gamma}_2{k}_{\mathrm{meteo}}\left(\left\{T,P,\omega, \mathrm{R}\mathrm{H}\right\}|h,\Big\{{T}^{\prime },{P}^{\prime },{\omega}^{\prime },\mathrm{R}{\mathrm{H}}^{\prime}\Big\}|{h}^{\prime}\right), $$

with respective variances $ {\gamma}_1 $ and $ {\gamma}_2 $ . We propose to parametrize the kernels using the Matérn family, which has been widely used to work with spatiotemporal data (Stein, Reference Stein1999). We denote $ {C}_{\nu } $ the Matérn- $ \nu $ covariance function with automatic relevance determination (that is independent lenghtscales for each entry). This choice of Matérn order $ \nu $ guarantees a desired level of regularity for the resulting GP. Details on the Matérn covariance functions are provided in Appendix A.

The spatiotemporal kernel is taken as the product kernel

(24) $$ {k}_{\mathrm{ST}}\left(\left\{t,\mathrm{la}\mathrm{t},\mathrm{lo}\mathrm{n}\right\},\left\{{t}^{\prime },\mathrm{la}{\mathrm{t}}^{\prime },\mathrm{lo}{\mathrm{n}}^{\prime}\right\}\right)={C}_{3/2}\left(|t-{t}^{\prime }|\right){C}_{3/2}\left(\mathrm{Hav}\left(\left[\begin{array}{c}\mathrm{la}\mathrm{t}\\ {}\mathrm{lo}\mathrm{n}\end{array}\right],\left[\begin{array}{c}\mathrm{la}{\mathrm{t}}^{\prime}\\ {}\mathrm{lo}{\mathrm{n}}^{\prime}\end{array}\right]\right)\right), $$

where $ \mathrm{Hav} $ denotes the Haversine distance—or great-circle distance—that determines the distance between two points on a sphere given their latitudes and longitudes. Each Matérn covariance admits a length scale, respectively denoted $ {\mathrm{\ell}}_t $ and $ {\mathrm{\ell}}_{\mathrm{lat},\mathrm{lon}} $ . Note that $ {C}_{3/2} $ may not be positive definite over Hav for certain lengthscales (Feragen and Hauberg, Reference Feragen and Hauberg2016; Borovitskiy et al., Reference Borovitskiy, Terenin, Mostowsky and Deisenroth2020), thus we set an upper bound on $ {\mathrm{\ell}}_{\mathrm{lat},\mathrm{lon}} $ to ensure positive-definiteness. The Matérn order $ \nu =3/2 $ ensures the GP is continuously differentiable with respect to time and space.

The spatiotemporal kernel (24) makes the predicted extinction smooth across time and space. Its product structure ensures that for distant times ( $ t,{t}^{\prime } $ ) or distance locations ((lat, lon), (lat $ ^{\prime } $ , lon $ ^{\prime } $ )), the spatiotemporal covariance vanishes to zero. Hence, the GP covariance of predictions distant in space or time will disregard spatiotemporal predictors. This prevents overfitting over spatiotemporal predictors.

The meteorological kernel is taken as an automatic relevance determination covariance

(25) $$ {k}_{\mathrm{meteo}}\left(\left\{T,P,\omega, \mathrm{RH}\right\}|h,\left\{{T}^{\prime },{P}^{\prime },{\omega}^{\prime },\mathrm{RH}^{\prime}\right\}|{h}^{\prime}\right)={C}_{1/2}\left(\left[\begin{array}{c}\mid T-{T}^{\prime}\mid \\ {}\mid P-{P}^{\prime}\mid \\ {}\mid \mathrm{RH}-\mathrm{RH}^{\prime}\mid \\ {}\mid \omega -{\omega}^{\prime}\mid \end{array}\right]\right), $$

where each dimension admits its own length scale parameter, respectively denoted $ {\mathrm{\ell}}_T,{\mathrm{\ell}}_P,{\mathrm{\ell}}_{\mathrm{RH}} $ , and $ {\mathrm{\ell}}_{\omega } $ (we permit ourselves to drop the conditioning on $ h $ and $ {h}^{\prime } $ for conciseness). The Matérn order $ \nu =1/2 $ ensures that the GP is continuous with respect to meteorological predictors.

The kernel on meteorological predictors (25) allows to GP to learn how meteorological predictors should modulate extinction. It introduces covariance between predicted extinctions if meteorological predictors are close, even if the spatiotemporal covariance is zero (because of the additive structure in (23)). In practice, we find that setting a unit variance for this kernel $ {\gamma}_2=1 $ helps prevent vanishing kernel signal when tuning the hyperparameters.

We include in Appendix B a comparison of this choice of covariance structure with a more conventional tensor form covariance and evaluates their respective predictive performance.

3.2.3. Sparse variational GP setup and tuning

We use $ p=200 $ inducing locations $ \mathbf{z} $ as an arbitrary choice (fewer or more inducing locations can be used following needs). They are initialized by randomly drawing samples at boundary layer altitude from $ \mathbf{x} $ . The variational distribution $ q\left(\mathbf{u}\right)=\mathcal{N}\left(\mathbf{u}|{\boldsymbol{\mu}}_{\mathbf{z}},{\boldsymbol{\Sigma}}_{\mathbf{z}}\right) $ is initialized with $ {\boldsymbol{\mu}}_{\mathbf{z}}=0 $ and $ {\boldsymbol{\Sigma}}_{\mathbf{z}}={\mathbf{I}}_p $ . We denote $ \Theta =\left\{{\boldsymbol{\mu}}_{\mathbf{z}},{\boldsymbol{\Sigma}}_{\mathbf{z}},\mathbf{z},\sigma, {\gamma}_1,{\gamma}_2,{\mathrm{\ell}}_t,{\mathrm{\ell}}_{\mathrm{lat},\mathrm{lon}},{\mathrm{\ell}}_T,{\mathrm{\ell}}_P,{\mathrm{\ell}}_{\mathrm{RH}},{\mathrm{\ell}}_{\omega}\right\} $ the joint set of variational and model hyperparameters. We maximize the ELBO objective (19) with respect to $ \Theta $ using the Adam optimizer (Kingma and Ba, Reference Kingma and Ba2015).

Since the dataset includes $ 8\times 96\times 192=147,456 $ atmospheric columns, evaluating the ELBO for all columns at each optimization step is computationally prohibitive. Instead, we use a stochastic gradient approach and sample random mini-batches of columns over which we compute an unbiased noisy estimate of the ELBO gradient. Whilst compromising statistical efficiency, the high variance estimates of gradients have regularizing properties that mitigate overfitting, and this approach is guaranteed to converge to a local optima for $ \Theta $ (Hoffman et al., Reference Hoffman, Blei, Wang and Paisley2013). To meet memory limitations, we set a mini-batch size of 64, and therefore choose a small learning rates of $ 0.01 $ to avoid taking large steps based on a few samples. When possible, a larger mini-batch size couple with a larger learning rate should allow to reach a better local optimum faster (Hoffman et al., Reference Hoffman, Blei, Wang and Paisley2013).

3.3.1. Exponential link

We choose an exponential link function $ \psi =\exp $ to warp the GP. This particular choice is mathematically convenient as it makes the prior $ \varphi \left(x|h\right) $ a log-normal random variable given by

(26) $$ \varphi \left(x|h\right)={e}^{f\left(x|h\right)-h/L}\sim \mathrm{\mathcal{L}}\mathcal{N}\left(m\left(x|h\right)-h/L,\sqrt{k\left(x|h,x|h\right)}\right). $$

A direct consequence is that we know the analytical expressions of its moments and quantiles, which alleviates the need for estimation procedures. This becomes of particular interest when making prediction, since the variational posterior that approximates $ \varphi \left(x|h\right)\mid \tau $ is also a log-normal distribution given by $ \mathrm{\mathcal{L}}\mathcal{N}\left({\overline{\mu}}_{x\mid h}-h/L,{\overline{\Sigma}}_{x\mid h}^{1/2}\right) $ following notations from Section 2.3.2.

Therefore, we obtain approximations of the posterior mean, variance, and quantiles in closed-form following

(27) $$ \unicode{x1D53C}\left[\varphi \left(x|h\right)|\tau \right]\approx \exp \left({\overline{\mu}}_{x\mid h}-\frac{h}{L}+\frac{1}{2}{\overline{\Sigma}}_{x\mid h}\right), $$

(28) $$ \mathrm{Var}\left(\varphi \left(x|h\right)|\tau \right)\approx \left({e}^{{\overline{\Sigma}}_{x\mid h}}-1\right)\exp \left(2{\overline{\mu}}_{x\mid h}-\frac{2h}{L}+{\overline{\Sigma}}_{x\mid h}\right), $$

(29) $$ {Q}_{\alpha}\left(\varphi \left(x|h\right)|\tau \right)\approx \exp \left({\overline{\mu}}_{x\mid h}-\frac{h}{L}+\sqrt{2{\overline{\Sigma}}_{x\mid h}}{\mathrm{\operatorname{erf}}}^{-1}\left(2\alpha -1\right)\right), $$

where $ {Q}_{\alpha } $ denotes the $ \alpha $ -quantile function and $ \operatorname{erf} $ denotes the error function. This allows for an efficient computation of the predicted mean extinction profile and confidence regions.

3.3.2. Rescaling the predicted profiles

The posterior distribution $ p\left(\mathbf{f}|\tau \right) $ updates the GP with the information that the mapping $ \varphi $ should (in expectation) vertically integrate to the observed AOD $ \tau $ . While instructive, this is a weak constraint on the prior and does not guarantee that the prediction will effectively integrate to the AOD of the atmospheric column.

To get stronger enforcement of column-integrated values and ensure predictions effectively integrate to the AOD in expectation, we propose to rescale the model against observed AOD. Because observations $ {\tau}_1,\dots, {\tau}_n $ are corrupted with high-frequency observational noise, we introduce a spatially smoothed version of the AOD observations that filters out this noise. By applying a spatial Gaussian smoothing filter to the ECHAM-HAM AOD, we obtain spatially smoothed versions of the observations which we denote $ {}^s{\tau}_1,\dots {,}^s{\tau}_n $ . A comparison of the fields is displayed in Figure 3.

Figure 3. Left: ECHAM-HAM 550 nm AOD. Right: Spatially smoothed ECHAM-HAM 550 nm AOD.

At inference time, we then substitute the predicted posterior of extinction with its rescaled version given by

(30) $$ {}^s\varphi \left({x}_i|h\right)\mid \tau =\frac{{}^s{\tau}_i}{\int_0^H\unicode{x1D53C}\left[\varphi \left({x}_i|h\right)|\tau \right]\hskip0.1em \mathrm{d}h}\varphi \left({x}_i|h\right)\mid \tau, $$

for the $ i $ th column. This ensures that, in expectation, $ {}^s\varphi \left({x}_i|h\right)\mid \tau $ effectively integrates along height to $ {}^s\tau $ . When choosing $ \psi =\exp $ , performing this rescale is equivalent to shifting the location of the posterior in (26) by $ {\log}^s{\tau}_i-\log {\int}_0^H\unicode{x1D53C}\left[\varphi \left({x}_i|h\right)|\tau \right]\hskip0.1em \mathrm{d}h $ . The integral is approximated with a trapezoidal scheme.

3.3.3. From $ {}^s\varphi \left(x|h\right)\mid \tau $ to $ {b}_{\mathrm{ext}}(h)\mid \tau $

The final stage of inference is to estimate the posterior distribution of $ {b}_{\mathrm{ext}} $ following the observation model (6). Because the chosen observation model conserves the mean, the posterior mean of $ {b}_{\mathrm{ext}}(h) $ is the same as the posterior mean of $ {}^s\varphi \left(x|h\right) $ . Higher-order moments however cannot be obtained analytically and will require estimation by sampling from $ {b}_{\mathrm{ext}}(h)\mid \tau $ . This can be achieved by sampling from $ {}^s\varphi \left(x|h\right)\mid \tau $ and then plugging the samples in the observation model (6) to sample from $ {\left.{b}_{\mathrm{ext}}(h)\right|}^s\varphi \left(x|h\right) $ .

The observation model scale parameter $ {\sigma}_{\mathrm{ext}} $ needs to be calibrated against observational noise. We set $ {\sigma}_{\mathrm{ext}} $ to the value that minimizes the integrated calibration index (ICI) given by

(31) $$ \mathrm{ICI}={\int}_0^1\mid {N}_{\sigma_{\mathrm{ext}}}\left(\alpha \right)-\left(1-\alpha \right)\mid \hskip0.1em \mathrm{d}\alpha, $$

where $ {N}_{\sigma_{\mathrm{ext}}}\left(\alpha \right) $ is the percentage of ECHAM-HAM extinction coefficient observations that fall within the $ \left(1-\alpha \right) $ credible interval of the distribution of $ {b}_{\mathrm{ext}}\mid \tau $ . The ICI is evaluated against $ {n}_{\mathrm{calib}}=200 $ random atmospheric columns from the dataset. These columns form a separate calibration set which is solely used to calibrate $ {\sigma}_{\mathrm{ext}} $ . This separate calibration set is negligible in size in comparison with the $ 147,456 $ atmospheric columns available in the simulation.

3.4.1. Baseline models

We use as a comparative baseline an idealized exponential profile, similar to the vertical profiles postulated in AOD retrieval algorithms. For a fair comparison, we rescale the baseline profiles such that their integrated values match the AOD, following the procedure described in Section 3.3.2. We complement this baseline with a log-normal observation model for the extinction coefficient and calibrate its scale parameter against a separate calibration following the procedure outlined in Section 3.3.3. This enables probabilistic evaluation of the extinction coefficient profiles predicted with the idealized baseline.

In addition, to examine the contribution to predictions of our modeling choices, we consider the following ablated versions of VAExtGP:

• • GP only: we remove the idealized exponential vertical contribution to the prior. The prior is specified as the exponential of a GP only, that is, $ \varphi \left(x|h\right)=\exp \left(f\left(x|h\right)\right) $ with $ f\sim \mathrm{GP}\left(m,k\right) $ .

• • ST only: we only use spatiotemporal covariates as input variables for the model, that is, $ x=\left(t,\mathrm{lat},\mathrm{lon}\right) $ . The kernel is set to $ k={k}_{\mathrm{ST}} $ .

• • Meteo only: we only use meteorological covariates as input variables for the model, that is, $ x=\left(T,P,\mathrm{RH},\omega \right) $ . The kernel is set to $ k={k}_{\mathrm{meteo}} $ .

3.4.2. Evaluation metrics

We use two categories of metrics to quantify the quality of the predicted vertical profiles: deterministic metrics, which compare only the posterior mean prediction to the groundtruth ECHAM-HAM extinction coefficient profiles, and probabilistic metrics which evaluate the entire posterior probability distribution against the extinction coefficients from ECHAM-HAM simulation and thus also assess the quality of the uncertainty quantification in the posteriors. The metrics used are detailed in Table 2.

Table 2. Evaluation metrics

Note. Deterministic metrics compare the predicted posterior mean $ \unicode{x1D53C}\left[{b}_{\mathrm{ext}}|\tau \right] $ to the ECHAM-HAM extinction coefficient; Probabilistic metrics evaluate the complete predicted posterior probability distribution of $ {b}_{\mathrm{ext}}\mid \tau $ against ECHAM-HAM extinction coefficient.

For each metric, we compute scores for pixels along the entire atmospheric columns and for pixels lying within the boundary layer only ( $ <2\hskip0.22em \mathrm{km} $ ). Scores are averaged across all pixels.