Aerosol optical depth disaggregation: toward global aerosol vertical profiles

Abstract

Aerosol-cloud interactions constitute the largest source of uncertainty in assessments of anthropogenic climate change. This uncertainty arises in part from the difficulty in measuring the vertical distributions of aerosols, and only sporadic vertically resolved observations are available. We often have to settle for less informative vertically aggregated proxies such as aerosol optical depth (AOD). In this work, we develop a framework for the vertical disaggregation of AOD into extinction profiles, that is, the measure of light extinction throughout an atmospheric column, using readily available vertically resolved meteorological predictors such as temperature, pressure, or relative humidity. Using Bayesian nonparametric modeling, we devise a simple Gaussian process prior over aerosol vertical profiles and update it with AOD observations to infer a distribution over vertical extinction profiles. To validate our approach, we use ECHAM-HAM aerosol-climate model data which offers self-consistent simulations of meteorological covariates, AOD, and extinction profiles. Our results show that, while very simple, our model is able to reconstruct realistic extinction profiles with well-calibrated uncertainty, outperforming by an order of magnitude the idealized baseline which is typically used in satellite AOD retrieval algorithms. In particular, the model demonstrates a faithful reconstruction of extinction patterns arising from aerosol water uptake in the boundary layer. Observations however suggest that other extinction patterns, due to aerosol mass concentration, particle size, and radiative properties, might be more challenging to capture and require additional vertically resolved predictors.

Impact Statement

Aerosol-cloud interactions (ACIs) represent the largest uncertainty in assessments of global warming, with uncertainty bounds that could offset global warming or double its effects. This uncertainty arises in part from the inability to observe aerosol amounts at an appropriate resolution. Instead aerosol optical depth (AOD)—observed globally through satellite products—is commonly used as a two-dimensional proxy for aerosols. Yet, obtaining a global estimate of aerosol vertical distribution in the atmosphere would be more informative and help better constrain ACIs uncertainty. Here, we show that using AOD and readily available vertically resolved meteorological predictors (temperature, pressure, relative humidity, updraft), a simple probabilistic model can already yield realistic vertical extinction profiles together with appropriate uncertainty quantification. This highlights that such simple modeling can benefit satellite products, leading to more accurate priors over aerosol vertical profiles.

1. Introduction

Aerosols are microscopic particles ( $ <5\hskip0.22em \unicode{x03BC} \mathrm{m} $ ) suspended in the atmosphere. They can come from natural sources (e.g., dust, sea salt) or be emitted by human activity (e.g., black carbon).

They influence the Earth’s energy budget with a negative radiative forcing that counteracts global warming from anthropogenic greenhouse gases emissions. A fraction of this negative forcing stems from aerosols’ direct scattering of incoming solar radiation (McCormick and Ludwig, 1967): this is the direct effect. A larger fraction of this forcing is due to their modulation of the radiative properties of clouds: this is the indirect effect. By acting as cloud condensation nuclei (CCN), additional aerosols can drive up the cloud droplet number while driving down the mean cloud droplet size. The resulting clouds are brighter, larger, and last longer (Twomey, 1977; Albrecht, 1989). They hence reflect more solar radiation and cool the Earth.

Unlike the direct effect which can, in principle, be well constrained (Watson-Parris et al., 2020), the magnitude of the forcing induced by the indirect effect is difficult to estimate. There are two reasons for this: (i) the physical processes underpinning aerosol-cloud interactions (ACIs) are not yet fully understood, which hinders the estimation of present-day forcing; (ii) the present-day forcing must be compared to the forcing at pre-industrial state, which is also particularly challenging to quantify (Carslaw et al., 2013). In fact, observational and model-based studies of ACIs still disagree on the magnitude of this forcing. As a result, ACIs contribute the largest uncertainty in present-day global warming (Masson-Delmotte et al., 2021).

To better estimate present-day forcing we require accurate, global measurements of CCN concentrations to assess radiative properties of clouds (Twomey, 1977; Albrecht, 1989). Unfortunately, measuring CCN concentrations can only be achieved in situ, and while field campaigns have already been undertaken to collect detailed CCN observations, these measurements are spatio-temporally sparse and provide insufficient constraints on the global distribution of aerosols (Andreae, 2009; Spracklen et al., 2011).

For lack of better observations, the aerosol optical depth (AOD) has been widely adopted as a first-order proxy of CCN concentration in an atmospheric column (Nakajima et al., 2001; Andreae, 2009; Clarke and Kapustin, 2010). The AOD is a measure of the extinction of solar radiation through an atmospheric column. It is denoted by $ \tau $ and defined at a given wavelength, time, latitude, and longitude by

(1) $$ \tau ={\int}_0^H{b}_{\mathrm{ext}}(h)\hskip0.1em \mathrm{d}h, $$

where $ {b}_{\mathrm{ext}} $ is the extinction coefficient¹ and the integral is taken over the height $ H $ of an atmospheric column. The AOD is appealing because it is routinely observed on a global scale by satellite products (Remer et al., 2005) which, unlike in situ observations, offer long-term global records.

However, the AOD is a column-integrated quantity and does not provide information on the vertical distribution of aerosols. This is limiting as their vertical distribution strongly influences both the magnitude and even the sign of the forcing induced by the indirect effect. For example, both modeling (Stier, 2016) and observational studies (Painemal et al., 2020) find AOD inadequate for assessing ACIs over vast subtropical ocean areas, which play a key role in determining the radiation balance of the Earth. Yet, in both studies, the vertically resolved aerosol extinction coefficient $ {b}_{\mathrm{ext}} $ shows a significantly higher correlation with CCN concentrations. Stier (2016) also highlights the importance of determining aerosol vertical distributions to provide stronger constraints on CCN at specific altitudes. In particular, the AOD fails to describe near-surface properties such as the concentration of aerosols in the boundary layer.

In this work, we propose to probe whether AOD observations can be used to constrain a global prior over aerosol vertical distributions. Formally, given an AOD observation $ \tau $ , we want to reconstruct the corresponding extinction coefficient profile $ {b}_{\mathrm{ext}} $ . This amounts to the task of reconstructing three-dimensional (3D) profiles using height-integrated two-dimensional (2D) observations and quantities that are easier to obtain in 3D, such as temperature and relative humidity.

Motivated by the study of cloud vertical structures, this task has been framed in the past as fully supervised learning (Leinonen et al., 2019), that is, assuming observations of groundtruth vertical profiles were available. Collecting vertically resolved observations of aerosols optical properties is also possible, using lidar-based remote sensing instruments (Winker et al., 2013) or ground-based sun-photometers (Holben et al., 1998). While valuable, these observations are however limited by their low spatiotemporal coverage and prone to corruption (e.g., low signal-to-noise ratio, clear-sky requirement). Compiling high-quality observational data of aerosol vertical profiles at a large scale is thus challenging, making fully supervised learning approaches inadequate.

Instead, we propose to draw from spatial disaggregation methodologies that only require observations at the aggregated level. Spatial disaggregation is the task of inferring subgrid details given coarse resolution spatial observations. Postulating an underlying fine-grained spatial field that aggregates into coarse observations, this problem can be framed as weakly supervised learning (Zhou, 2017) with aggregated targets. While existing works (Law et al., 2018; Tanaka et al., 2019; Yousefi et al., 2019; Tanskanen and Longi, 2020; Zhang et al., 2020) have only considered aggregation processes happening on a 2D field, this rationale can be extended to disaggregate quantities along a third dimension—height. Since $ \tau $ corresponds to the vertical integration of $ {b}_{\mathrm{ext}} $ , we propose to frame the reconstruction of aerosol vertical profiles as the vertical disaggregation of AOD observations.

Using Gaussian processes (GPs) (Rasmussen and Williams, 2005), we design a Bayesian model that maps vertically resolved meteorological variables (e.g., pressure, relative humidity) to a probabilistic estimate of the extinction coefficient that integrates into the AOD. The model formulation is simple and makes assumptions explicit, hence granting control and interpretability over predictions while offering built-in uncertainty quantification.

In order to be able to fully validate the proposed methodology, we use ECHAM-HAM global aerosol-climate model simulation data (Stier et al., 2005, 2007; Zhang et al., 2012). While our prime motivation is to reconstruct $ {b}_{\mathrm{ext}} $ from satellite observations of AOD, the intricacies of combining measurements from different instruments make it challenging to validate any proposed methodology. On the other hand, ECHAM-HAM is a self-consistent climate model that offers readily available aerosol vertical profiles, and is better suited for model development. We demonstrate our model is able to reconstruct natural patterns that arise in aerosol vertical distribution, in particular in the boundary layer. We show that very simple and readily available meteorological predictors suffice to obtain a good estimation of the extinction coefficient.

The aims of this study are as follows:

• • Formulate a probabilistic vertical disaggregation methodology for the task of reconstructing aerosol vertical profiles using readily available vertically resolved covariates.

• • Validate the proposed methodology using climate model data—where access to groundtruth extinction coefficient profiles enables quantitative evaluation.

Section 2 outlines the design of the vertical disaggregation model and describes the inference procedure as well as the hyperparameter selection. Section 3 describes the dataset and experimental setup used to validate the model. Section 4 presents the experimental results and Section 5 discusses, while introducing avenues of future research.

2. Model design

In this section, we first outline the design of the prior distribution we place on the extinction coefficient profile using GPs. We then describe the observation model that connects our prior to the observations of the AOD. Finally, we present how the posterior inference over the extinction profiles is conducted.

2.1. Design of the prior

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., 2017). 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., 2017; Li et al., 2020). 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., 2004; Small et al., 2009; Chen et al., 2014; Douglas and L’Ecuyer, 2020)—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, 1936; Twomey, 1974): 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., 2017; Jesson et al., 2022). 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., 2021), 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., 2022). 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, 2005) 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., 2016). 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, 1999; Rasmussen and Williams, 2005).

2.2. Design of the observation models

In this section, we establish a connection between the extinction coefficient prior $ \varphi \left(x|h\right) $ (constructed in Section 2.1) and the observed values of extinction coefficient and AOD. To simplify notations, we focus on a single air column with a height of $ H $ , where we have observations of the extinction coefficient $ {b}_{\mathrm{ext}}(h) $ and the AOD $ \tau $ .

Ideally, we would like to achieve a perfect match between $ {b}_{\mathrm{ext}}(h) $ and $ \varphi \left(x|h\right) $ , as well as between $ \tau $ and $ {\int}_0^H\varphi \left(x|h\right) $ . Unfortunately, this is unrealistic because observations are prone to measurement error, which must critically be accounted for. To address this, we formulate probabilistic observation models for $ {b}_{\mathrm{ext}} $ and $ \tau $ in what follows.

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., 2000). 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., 1998), 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 mean² 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., 2005) or assumptions made in the AOD retrieval algorithm (Mielonen et al., 2011; Wu et al., 2016). 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. Tuning and posterior inference

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, 2009; Matthews et al., 2016; Leibfried et al., 2020) to substitute this intractable inference problem with a tractable optimization problem. In addition, we make the variational approximation sparse (Titsias, 2009) 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., 2015). 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., 2019). As it is standard in sparse variational GPs (Titsias, 2009), we also learn the inducing locations $ \mathbf{z} $ .

3. Experimental setup

Our motivating application is to reconstruct aerosol vertical profiles from remote-sensing AOD observations. However, to validate the proposed methodology, we also need to observe extinction coefficient profiles. While the latter can be collected from vertically resolved remote sensing instruments (Holben et al., 1998; Winker et al., 2013), the collocation of measurements from different devices raises non-trivial questions regarding the self-consistency of the collocated dataset—and thus the validation procedure. In contrast, climate models are self-consistent and offer readily available simulations of gridded AOD and extinction coefficient profiles. We hence propose to evaluate our model using observations from a global climate model, the ECHAM-HAM global aerosol-climate model simulation (Stier et al., 2005, 2007; Zhang et al., 2012).

The objective of this experiment is to apply our model to the vertical disaggregation of the AOD simulated with ECHAM-HAM, and compare our prediction to the extinction coefficient simulated with ECHAM-HAM. In this section, we first describe the ECHAM-HAM dataset used in the experiments. We then outline the model setup and the inference procedure. Finally, we present the evaluation metrics used to quantify the quality of the reconstructed vertical profiles.

3.1. Dataset

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., 2005). 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, 2022); (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, 2022). 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., 2005 and 532 nm for the CALIOP extinction; Winker et al., 2013).

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. Model setup and tuning

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., 2007; Kahn and Gaitley, 2015; Laszlo and Liu, 2016). 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, 1999). 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, 2016; Borovitskiy et al., 2020), 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, 2015).

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., 2013). 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., 2013).

3.3. Inference procedure

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. Evaluation procedure

We compare the predicted extinction coefficient profile $ {b}_{\mathrm{ext}}\mid \tau $ to the extinction coefficient from ECHAM-HAM simulations. For clarity, we refer to our proposed methodology as VAExtGP (for Variational Aerosol Extinction GP) in the presentation of experimental results.

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.

4. Results

4.1. Comparison to an idealized exponential baseline

Table 3 shows that the posterior mean profile predicted with the proposed method outperforms the idealized exponential baseline across deterministic metrics. We note that evaluating over the entire column consistently yields better scores compared to evaluating only over the boundary layer. This is to be expected as the extinction coefficient outside the boundary layer tends to vanish to zero and most of the variability occurs within the boundary layer.

Table 3. Scores of the idealized exponential baseline and VAExtGP for the task of predicting ECHAM-HAM extinction profiles

Note. “Entire column” means scores are computed for every altitude level; “Boundary layer” means scores are computed for altitude levels of the boundary level only ( $ <2\hskip0.22em \mathrm{km} $ ); the best scores for each metric and region are highlighted in bold; we report 1 standard deviation.

^a For Idealized, we report the marginal log-likelihood of ECHAM-HAM extinction under the predicted profiles instead of the ELBO. This is an upper bound to the ELBO, which could possibly be lower.

The improvement in RMSE and MAE against the idealized baseline is slightly more significant when computed within the boundary layer rather than over the entire column. This supports the idea that VAExtGP improves upon an idealized exponential profile at predicting extinction since the boundary layer concentrates most extinction variability. The improvement is greatest in Bias and Bias98, where the proposed model consistently improves over the idealized baseline by one or multiple orders of magnitude. In particular, the improvement in Bias98 suggests that VAExtGP better captures extreme extinction values.

When considering the evaluation over the entire column, the ELBO and ICI scores are comparable for both methods. However, within the boundary layer, VAExtGP outperforms the baseline, by a substantial margin in terms of ICI. This indicates that, within the boundary layer, the predicted posterior probability distribution of VAExtGP better captures the ECHAM-HAM extinction variability. The 95% calibration scores are similar for both methods, with a slight advantage for the idealized baseline.

4.2. Reconstructed vertical profiles

Figure 4 displays slices for a fixed latitude of the vertically resolved predictors used, the ECHAM-HAM extinction coefficient, and the predicted extinction coefficient with VAExtGP. For comparison, the prediction at the same latitude with the idealized exponential baseline and additional predictions with VAExtGP for different latitudes are provided in Appendix C.

Figure 4. Vertical slices at latitude 51.29° of meteorological predictors ( $ T,P,\mathrm{RH},\omega $ ), groundtruth extinction coefficient, predicted extinction coefficient posterior mean, 2.5% and 97.5% quantiles of the predicted extinction coefficient posterior distribution.

We observe that our predicted mean profile is able to reconstruct extinction patterns that are visually very similar to the groundtruth extinction coefficient. In comparison with the profiled predicted with the idealized exponential baseline in Figure C1, the extinction profiles predicted with our method look much more realistic. This is encouraging considering that the only aerosol optical property used is the AOD. We also observe that extreme extinction coefficient values appear to be well captured within the 95% confidence region of the predicted posterior distribution.

Since aerosol water uptake is related to relative humidity, we observe a good capacity to recover patterns corresponding to extinction due to aerosol swelling at low altitudes in the boundary layer. This is particularly visible over the Pacific (longitudes $ {144}^{\circ }-{228}^{\circ } $ ) and the Atlantic (longitudes $ {305}^{\circ }-{350}^{\circ } $ ), where the predicted patterns of extinction tend to best reproduce the ECHAM-HAM extinction patterns.

However, the predicted mean profile faces greater difficulties to accurately reproduce the extinction patterns for longitudes $ {0}^{\circ }-{100}^{\circ } $ . This is particularly noticeable in Figure 4 around longitude $ {71.4}^{\circ } $ nearby Astana in Kazakhstan, and around longitude $ {17.3}^{\circ } $ in the Silesia region in Poland. The geographical locations suggest these are extinction patterns induced by human activity, which cannot be captured solely by meteorological variables. We conjecture these patterns are imputable to aerosols mass concentration, particles size, and radiative properties. These properties, unlike extinction due to swelling, are more complex and cannot be fully characterized by relative humidity, temperature, pressure, and updraft. It is worth highlighting that, whilst difficult to predict, these extinctions pattern appear nonetheless to be well captured within the 95% confidence region of the predicted posterior distribution.

In general, the posterior mean predictions tend to be smoother than the groundtruth extinction coefficient. The smoothing effect is a regularizing property of the GP prior, resulting in an overestimation of the extinction coefficient in certain regions with low extinction. For example, around longitude $ {200}^{\circ } $ in the Pacific, the ECHAM-HAM extinction in Figure 4 displays an extinction pocket with a sharp limit around altitude $ 1\hskip0.22em \mathrm{km} $ , above which extinction is virtually absent. In contrast, our predicted posterior mean is more diffuse and tends to spread out lightly above $ 1\hskip0.22em \mathrm{km} $ .

The density plots in Figure 5 support this observation, as a large mass of extinction coefficient around $ {10}^{-7}\hskip0.1em {\mathrm{m}}^{-1} $ tends to be overestimated by an order of magnitude. However, this trend is significantly reduced when focusing on the boundary layer, where most of the mass lies around the axis $ y=x $ . This suggests that most of the overestimation occurs above the boundary layer for low extinction ( $ <{10}^{-6}\hskip0.1em {\mathrm{m}}^{-1} $ ). Within the boundary layer, however, the mean predictions are reasonably aligned with the groundtruth for high extinction coefficient values ( $ >{10}^{-5}\hskip0.1em {\mathrm{m}}^{-1} $ ). Although low extinction coefficient values are less dominant in the boundary layer, they also tend to be overestimated. We attribute this to the smoothing effect of the GP prior, which is prone to overestimation between extinction pockets.

Figure 5. Density plots of groundtruth extinction coefficient values against predicted posterior mean extinction coefficient; Left: plotted for the entire column; Right: plotted for the boundary layer only; density plots are computed on a random subset on a random subset of 1000 samples drawn for the entire column (left) and in the boundary layer (right).

4.3. Ablation study

We examine the contribution to predictions of our modeling choices by comparing the scores obtained with VAExtGP to the scores obtained with the three ablated models described in Section 3.4.1. Figure 6 shows the relative changes in scores with the ablated models, defined as

(32) $$ \%\mathrm{Relative}\ \mathrm{change}=\frac{\mathrm{Score}\; VAExtGP-\mathrm{Score}\ \mathrm{ablated}\ \mathrm{model}}{\mathrm{Score}\; VAExtGP}\times 100, $$

where a positive value (blue) corresponds to an improvement for the ablated model, and a negative value (red) corresponds to a decrease in performance for the ablated model. Transformations of Corr, Bias, Bias98, ELBO, and Calib95 are taken such that smaller values correspond to better performance, and are detailed in Appendix C.

Figure 6. Relative changes in scores with respect to VAExtGP for the 3 ablated models: “GP only” (no idealized exponential term $ {e}^{-h/L} $ in the prior), “ST only” (only spatiotemporal covariates in the input variable) and “Meteo only” (only meteorological covariates in the input variable). Red/Blue indicates that the performance is on average worse/better for the ablated model. We report 1 standard deviation.

Overall, we find that excluding a component from the model leads to a decrease in performance across most metrics. Specifically, when using only a GP without the idealized exponential term in the prior (GP only), only spatiotemporal covariates (ST only), or only meteorological covariates (Meteo only), the model’s performance is adversely affected.

For deterministic metrics, the ablation of meteorological covariates (ST only) has a particularly pronounced impact, with drastic drop in Bias and Bias98. This highlights the key role of learning the mapping from meteorological variables to extinction profiles. For probabilistic metrics, we observe that the ablation of spatiotemporal covariates (Meteo only) leads to the largest decrease in scores, indicating that incorporating spatiotemporal information plays a significant role in calibrating the predicted posterior distribution. Furthermore, we note that whilst the performance of the GP only model is generally inferior to that of VAExt, emphasizing the importance of encoding prior knowledge about aerosol vertical distribution through the idealized exponential profile, it is also the ablated model that least degrades performance. This highlights the remarkable adaptivity of GPs, and supports their use to learn non-trivial relationship in data.

In conclusion, these results demonstrate the significance of incorporating both spatiotemporal and meteorological covariates, as well as the idealized exponential profile, in our proposed model.

4.4. Feature importance analysis

In this section, we investigate the role predictors in VAExtGP play in the prediction of the extinction profiles. We begin by examining the model hyperparameters to gain an understanding of how do predictors overall contribute to the covariance structure. We then ask the following questions:

• • How do meteorological predictors modulate the predicted extinction?

• • Do temperature and pressure provide valuable meteorological information for prediction or are they just proxies for height?

To address these questions, we propose a more detailed feature importance analysis that elucidates the contribution of predictors to individual predictions: Shapley values (SVs) (Shapley, 1953). Originally introduced as a game-theoretic concept, SVs have been widely adopted in machine learning to design local feature importance explanation models (Štrumbelj and Kononenko, 2014; Lundberg and Lee, 2017; Lundberg et al., 2018; Ghorbani and Zou, 2019). With SVs, we can estimate for each individual prediction, how much each predictor contributed, and hence explain predictions. In our study, we use the KernelSHAP SV-based explanation model (Lundberg and Lee, 2017).

4.4.1. Global analysis with model hyperparameters

Table 4 reports the hyperparameters values after tuning them to maximize the ELBO objective as described in Section 2.3.3. These values provide a global description of the importance of predictors contribution to the covariance structure. $ {\gamma}_1 $ and $ {\gamma}_2 $ respectively correspond to the variances of the spatiotemporal and meteorological kernels from (23). A larger variance hyperparameter indicates a greater contribution to the covariance. The lengthscales $ {\mathrm{\ell}}_{\mathrm{lat},\mathrm{lon}},{\mathrm{\ell}}_T,{\mathrm{\ell}}_P,{\mathrm{\ell}}_{\mathrm{RH}} $ , and $ {\mathrm{\ell}}_{\omega } $ correspond to the lengthscales of each standardized predictor in the spatiotemporal kernel (24) and (25). A smaller length scale parameter means that the covariance is more sensitive to small variations in the predictors.

Table 4. Hyperparameter values for VAExtGP tuned for standardized input variables following the procedure described in Section 2.3.3

Note. Length scale is unitless and apply to standardized predictors; hyperparameters are tuned using 10 different initialization seeds; we report 1 standard deviation; $ \dagger $ indicates the value is fixed in the model and not tuned against data.

The value of $ {\gamma}_1 $ indicates that meteorological predictors globally contribute more to the covariance than spatiotemporal predictors. The large value of the temporal length scale $ {\mathrm{\ell}}_t $ (greater than 5) suggests a minimal sensitivity of the model to time. This is certainly because the dataset covers a single day and therefore exhibits limited temporal variability. In contrast, the small value of the spatial length scale $ {\mathrm{\ell}}_{\mathrm{lat},\mathrm{lon}} $ indicates a high sensitivity to spatial location and suggests that spatial correlation occurs within a restricted spatial range. This likely arises from the spatial heterogeneity of aerosol concentration, resulting in localized, short-range spatial patterns.

Regarding meteorological predictors, the length scales indicate that—in general—predictions are more sensitive to variations in temperatures and less sensitive to variations in updraft. In fact, the large value of the updraft length scale $ {\mathrm{\ell}}_{\omega } $ indicates, in general, a minimal correlation decay over the range of inputs.

4.4.2. How do meteorological predictors modulate the predicted extinction?

To gain a refined understanding of how meteorological predictors affect predictions, we conduct a SV analysis of their contribution to individual predictions. Since the meteorological predictors only modulate predictions through the weighting function, we only study the impact meteorological predictors have on the weight $ w\left(x|h\right) $ predicted for the posterior mean response. We also focus on samples within the boundary layer since this is where most of the variability happens. Figure 7 displays the mean absolute SV obtained for each meteorological predictor and the contribution each predictor has on a randomly selected set of individual predictions.

Figure 7. Left: Mean absolute Shapley values of meteorological predictors in the boundary layer; Right: Beeswarm plot of Shapley values of individual predictions in the boundary layer for each meteorological predictor; red/blue dots indicate a sample where the predictor value lies on the upper/lower end of its distribution; Shapley values are computed on a random subset of 2000 samples in the boundary layer.

We observe that in absolute mean, temperature is the factor that most influences predictions of the weighting function in the boundary layer. This is sensible as temperature is a smooth spatial field that gradually decreases with altitude, making $ T $ characteristic of the altitude level. Since the altitude typically correlates with extinction values, $ T $ becomes an informative proxy of extinction. This is corroborated by the plot of contributions to individual predictions: greater temperature (red)—lower altitudes—contribute to a greater weight $ w\left(x|h\right) $ and hence greater extinction prediction while lower temperature (blue)—higher altitudes—yield lower impact on the weight and hence lower extinction prediction. It is interesting to notice that while pressure is also a proxy of altitude just like temperature, its mean absolute contribution is lower and the individual contribution plot does not display a clear identifiable trend. This suggests that the information conveyed by pressure might be redundant with the information conveyed by temperature, and could possibly be discarded from the model without loss of predictive performance.

Whilst the mean absolute contribution of relative humidity seems marginal compared to temperature, we notice in individual contributions that when RH contributes to a greater weight $ w\left(x|h\right) $ (points lying on the right), it consistently corresponds to locations with high RH. This supports the observation that our model captures extinction due to aerosol water uptake. Similarly, the mean absolute contribution of $ \omega $ is relatively small, but a clear trend appears in the individual contributions: a positive updraft consistently corresponds to a positive contribution of $ \omega $ to the weight function, whereas a negative updraft corresponds to a negative contribution to the weight function. This is consistent with intuition as strong updrafts lead to the humidification of aerosols, and an increase in water content in aerosols increases in turn the extinction due to aerosol water uptake. Conversely, a negative updraft should lead to lower humidification of aerosols, and therefore a decrease in extinction due to aerosols water uptake.

4.4.3. Temperature and pressure: proxies for height?

Temperature and pressure carry significant information about height, and may therefore inform VAExtGP primarily as a proxy for height. Indeed, the vertical slices of $ T $ and $ P $ in Figure 4 exhibit an approximately linear behavior with height. However, for a fixed height, temperature and pressure also exhibit variations along longitude, and may therefore provide additional information about meteorological conditions, beyond serving as a proxy for height. To investigate this, we conduct a SV analysis of the contribution of height, temperature, and pressure to the predicted logarithm of the posterior mean $ \log \unicode{x1D53C}\left[\varphi \left(x|h\right)|\tau \right] $ . Figure 8 shows the reported SVs for the same vertical slice as Figure 4.

Figure 8. Vertical slices at latitude 51.29° (same as Figure 4) of the Shapley values for height h, temperature T, and pressure P; the Shapley values are computed for predicting the logarithm of the posterior mean, that is, $ \log \unicode{x1D53C}\left[\varphi \left(x|h\right)|\tau \right] $ ; green/pink indicates a sample where the predictor has contributed to a lower/higher predicted extinction.

The contribution of $ h $ to predicted extinction is clearly visible, where greater altitude contributes to a lower predicted extinction and lower altitude contributes to a greater extinction. This corresponds to the domain knowledge we encode within the prior by using an idealized exponential component. The contribution of $ T $ displays a similar mode, in agreement with the analysis from Section 4.4.2: at high altitude, lower temperatures contribute to lower predicted extinction, whereas at low altitudes, higher temperatures contribute to greater predicted extinction. The contribution of $ P $ is in general marginal in comparison, which aligns with the conjecture that $ P $ carries redundant information and may be a confounding feature which can be discarded.

In the stratosphere (altitude $ >10\;\mathrm{km} $ ), the contribution of $ T $ becomes marginal in comparison with the contribution of $ h $ . This suggests that the contribution of the exponential profile dominates at high altitudes. In the troposphere (altitude $ <10\hskip0.24em \mathrm{km} $ ), the contribution of $ T $ dominates, indicating that temperature plays a key role in predictions at lower altitudes. This is particularly visible in the boundary layer (altitude $ <2-3\hskip0.24em \mathrm{km} $ ). Importantly, the temperature SVs display spatial heterogeneity within the boundary layer, where the contribution at a fixed altitude varies across longitude. This is different from the height SVs, and suggests that $ T $ informs the model with valuable meteorological information, beyond being a proxy for height.

To evaluate the contribution of $ T $ and $ P $ to prediction performance, we conduct ablation experiments where we train and evaluate VAExtGP by (i) removing $ P $ and (ii) removing both $ T,P $ from the input meteorological predictors. Table 5 shows the percentages of relative change in score in comparison with the entire VAExtGP model. We follow the same conventions than in Section 4.3 and take transformations of Corr, Bias, Bias98, ELBO, and Calib95 such that smaller values correspond to better performance.

Table 5. Percentages of relative change in scores in with respect to VAExtGP when removing $ P $ and $ \left(T,P\right) $ from the input meteorological predictors

Note. Red/Blue indicates that the performance is worse/better for the ablated model with statistical significance at $ p<0.05 $ level. We report 1 standard deviation.

We find that: (i) predictive performance is overall improved by removing $ P $ . This supports the conjecture that $ P $ is a confounding factor, which carries redundant information for the prediction of extinction, and thereby is a source of noise in predictions. (ii) Removing both $ T $ and $ P $ overall harms predictive performance. This supports the observation that, whilst strongly correlated with height, temperature and pressure provide additional predictive utility.

5. Conclusion and outlooks

In this work, we introduce a GP-based methodology to vertically disaggregate the AOD using simple vertically resolved meteorological predictors such as temperature, pressure, or relative humidity. Our approach emphasizes uncertainty quantification using a Bayesian formalism. A successful application of our methodology to the vertical disaggregation of ECHAM-HAM simulated AOD is demonstrated. Our model outperforms an idealized baseline and displays capacity to recover realistic extinction patterns, in particular for extinction patterns arising from aerosol swelling in the boundary layer.

The simplicity of the explicit modeling assumptions grants better control and interpretability over the model and makes the choice of analysis strategy less subjective. While such simplicity can never account for the complex phenomena underpinning ACIs, this is balanced by a particular emphasis on the quantification of epistemic uncertainty through a principled Bayesian formalism.

Naturally, the modeling assumptions can be adapted to reflect different modeling choices. For example, where more relevant vertically resolved quantities can be obtained, they can seamlessly be incorporated into the predictors. The idealized exponential component of the prior can be replaced by an idealized Gaussian profile to reflect the assumptions made by different AOD retrieval algorithms (Wu et al., 2017; Li et al., 2020). The kernel design is also flexible and allows the user to specify a covariance structure that best accounts for pairwise dependencies of the predictors. Similarly, the positive link function, integration scheme, or number of inducing locations can be modified to fit specific needs.

Experiments demonstrate the proposed model is able to realistically reproduce vertical structures of extinction. In particular—for the chosen set of meteorological predictors—we are able to reliably predict extinction due to water vapor in the boundary layer. Because of the smoothness bestowed by the GP prior, the model tends to overestimate low extinction above the boundary layer and between extinction packets in the boundary layer. It is conjectured that the remaining unexplained extinctions patterns can be attributed to aerosols mass concentration, particles size, and radiative properties, which are more challenging to model and would require additional vertically resolved covariates. This constitutes an important area for future work.

Regarding methodology, several extensions remain to be explored. For example, while we only consider aerosol extinction at a single wavelength, the AOD is in general retrieved for multiple wavelengths ranging from $ 470\hskip0.22em \mathrm{nm} $ to $ 870\hskip0.22em \mathrm{nm} $ (Remer et al., 2005). By making the GP a function of the wavelength $ f\left(\lambda, x|h\right) $ , it should be possible to leverage multiple AOD observations by introducing a notion of functional smoothness across the electromagnetic spectrum. Another exciting direction would be to allow the use of AOD observations which are not matched with vertically resolved predictors. This can be achieved by adapting the work of Chau et al. (2021) to our model, using a globally observed 2D field that would mediate the learning between unmatched AOD and vertically resolved predictors. Such addition would be particularly welcome as it would allow prediction of aerosol vertical profiles even at locations where AOD is not observed but vertically resolved predictors are.

The proposed method can also be useful to diagnose the extinction profiles produced by Earth system models (ESMs). For example, by training and validating our model on ESMs, one can produce a predictive model that learns ESM extinction characteristics. Predictions from this model can then be evaluated against real extinction observations, and thereby provide a quantitative indicator of whether ESM extinction characteristics match the characteristics of observed extinction.

Finally, a different exciting use case of the proposed methodology can be considered. If working with climate model data only, then one can choose to use any vertically resolved aerosol tracers as predictors. These include mass and number concentrations which are simulated for several aerosol modes (nucleation, aitken, coarse, and accumulation) and species (dust, sea salt, sulfate, black, and organic carbon). Leveraging the interpretability of the GP covariance function, the kernel hyperparameters can then indicate which mode/specie did contribute to predictions, hence providing a tool to better understand factors that modulate the optical properties of aerosols.

In future work, we intend to apply our model to AOD arising from satellite observations and meteorological predictors from reanalysis data. Since observations of groundtruth extinction coefficients are not available in 2D satellite products, we intend to collocate observations from MODIS 2D AOD products (Remer et al., 2005) with CALIOP vertical lidar measurements (Winker et al., 2013) to validate the model.