the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Introduction of the DISAMAR radiative transfer model: determining instrument specifications and analysing methods for atmospheric retrieval (version 4.1.5)
Johan F. de Haan
Ping Wang
Maarten Sneep
J. Pepijn Veefkind
Piet Stammes
DISAMAR (determining instrument specifications and analysing methods for atmospheric retrieval) is a computer model developed to simulate retrievals of properties of atmospheric trace gases, aerosols, clouds, and the ground surface from passive remote sensing observations in a wavelength range from 270 to 2400 nm. It is being used for the TROPOMI/Sentinel5P and Sentinel4/5 missions to derive Level1b product specifications. DISAMAR uses the doubling–adding method and the layerbased orders of scattering method for radiative transfer calculations. It can perform retrievals using three different approaches: optimal estimation (OE), differential optical absorption spectroscopy (DOAS), and the combination of DOAS and OE, called DISMAS (differential and smooth absorption separated). The derivatives, which are needed in the OE and DISMAS retrievals, are derived in a semianalytical way from the adding formulae. DISAMAR uses planeparallel homogeneous atmospheric layers with a pseudospherical correction for large solar zenith angles. DISAMAR has various novel features and diverse retrieval possibilities, such as retrieving aerosol layer heights and ozone vertical profiles. This paper provides an overview of the DISAMAR model version 4.1.5 without treating all the details. We focus on the principle of the layerbased orders of scattering method, the calculation of the semianalytical derivatives, and the DISMAS retrieval method, and it is to our knowledge the first time that these methods are described. We demonstrate some applications of DISMAS and the derivatives.
DISAMAR stands for determining instrument specifications and analysing methods for atmospheric retrieval. It is a computer code written in Fortran 90 to simulate retrievals of atmospheric components like trace gases, aerosols, and clouds and properties of the ground surface from passive remote sensing observations of the Earth. The wavelength range considered is from 270 to 2400 nm. The development of DISAMAR started around the year 2000 but made use of the heritage of doubling–adding code from the 1980s (De Haan et al., 1987). Initially, DISAMAR was developed to derive Level1b radiance requirements given specific Level2 geophysical product requirements for satellite instruments. In order to derive Level1b requirements for a given requirement for a Level2 product, the model has to simulate instrument features and perform retrievals for the Level2 product with proper error estimation. From the European satellite missions Sentinel5P (with TROPOMI on board Veefkind et al., 2012), Sentinel4, and Sentinel5, atmospheric trace gases, cloud, and aerosol products have to be generated from the measured radiance spectra. In principle, it is possible to simulate Level1b spectra using a radiative transfer (RT) model to modify the simulated spectra with instrument features, like noise, and then run a retrieval algorithm. However, this can introduce inconsistencies between the RT model simulation and the retrieval algorithm, which makes the error estimate of the Level2 product less accurate. With DISAMAR providing the endtoend calculation for all the steps from instrument response to retrieved geophysical product, a high degree of internal consistency is achieved. For example, DISAMAR has been used to calculate Level1b requirements for preparation of the Sentinel4, Sentinel5, Sentinel5P and other missions (e.g. De Haan, 2009, 2010, 2012; Sanders et al., 2015; Nanda et al., 2018).
DISAMAR contains a radiative transfer module to simulate radiance, irradiance, and sunnormalized radiance (or reflectance). The radiative transfer is based on the doubling–adding method and includes a more efficient variant, called layerbased orders of scattering (LABOS). Polarization (Stokes fourvectors and 4×4 Mueller matrices) and rotational Raman scattering (RRS) are included in DISAMAR. The implementation of RRS is based on the publications by Chance and Spurr (1997) and Stam et al. (2002). The derivatives of the radiance or reflectance with respect to the retrieved parameters (also called weighting functions or Jacobians) that are needed for the retrievals are calculated using efficient semianalytical expressions. The atmosphere is assumed to be planeparallel with an option for pseudospherical correction. Three types of retrievals are available, namely optimal estimation (OE) (Rodgers, 2000), DOAS (differential optical absorption spectroscopy) (Platt, 1994), and DISMAS (differential and smooth absorption separated), a version of OE based on a DOASlike approach. The choice of DOAS and OE retrievals is driven by the Level2 retrieval algorithms for atmospheric trace gas columns (using DOAS) and ozone profiles (using OE).
A cloud or aerosol layer can be modelled by a simple reflecting Lambertian surface or by a layer of scattering particles with a Henyey–Greenstein phase function or a phase function using expansion coefficients that are read from file (e.g. Mie scattering particles). The surface below the atmosphere is a Lambertian reflector with an albedo that can vary with wavelength. In addition, surface emissions can be included to simulate suninduced fluorescence (Frankenberg et al., 2011).
The applications of DISAMAR are retrievals of the ozone profile from the ultraviolet (UV), the total column of NO_{2} from the visible (VIS), cloud properties (height, fraction, optical thickness) from the O_{2} Aband in the nearinfrared (NIR), or CO_{2} from the shortwave infrared (SWIR). The simulated measurements can be modified by adding noise, offsets to the wavelength scales, and offsets to the radiance and solar irradiance. By evaluating the effects of these modifications on retrieved parameters, requirements on the Level1b measured spectra can be derived. By changing the a priori error of auxiliary data (e.g. surface albedo) the required accuracy of auxiliary data can be determined.
Due to the timeconsuming linebyline calculations, DISAMAR is not suitable for fast computation or applications in numerical weather prediction (NWP) models. We would recommend RTTOV (Saunders et al., 2018) and CRTM (e.g. Lu et al., 2021; Karpowicz et al., 2022; Stegmann et al., 2022) as fast radiative transfer models for NWP applications. DISAMAR has actually been used to benchmark RTTOV in the UV–visible wavelength range.
The advantage of DISAMAR is that it includes diverse instrument features, simulations, and retrieval algorithms in one code. For such a complex model, it is not possible to cover all details in one paper. We therefore highlight the novel and unique features of DISAMAR (version 4.1.5), which are listed in Sect. 2. Section 3 presents the forward simulations. The retrieval part is described in Sect. 4. Some applications are shown in Sect. 5 as examples. Section 6 presents a discussion and conclusions.
DISAMAR uses a separate altitude grid for the radiative transfer calculations, which is independent of the grid used for specifying the atmospheric properties. This makes it possible to deal with strong vertical gradients in the radiation field, e.g. near the top of clouds.
The altitude grid, wavelength grid, and polar angles are all defined using Gauss–Legendre division points. This makes the integration over altitude, wavelength, or polar angle more accurate than the integration with an equidistant grid using a similar number of grid points.
The number of streams used for integration over polar angles when multiple scattering is involved can be arbitrarily large. The number of coefficients used for the expansion of the phase function in Legendre functions can also be arbitrarily large. All of these features make the DISAMAR calculations more accurate.
DISAMAR provides not just the radiance spectrum of the backscattered sunlight but also the derivatives with respect to the elements of the state vector. These derivatives are essential in OE and are used to find the solution in an iterative manner. They are also used to determine the error covariance matrix, gain vectors, and the averaging kernel (Rodgers, 2000). In DISAMAR all derivatives are calculated in a semianalytical manner, although algorithmic differentiation can be used to evaluate the derivatives (Griewank and Walther, 2008). The calculation of the derivatives in DISAMAR takes no more than 10 %–30 % of the total calculation time.
Because DISAMAR was originally developed to determine instrument specifications and analyse methods for atmospheric retrieval, special emphasis is given to the error information of the retrieved products and instrumental features. The main retrieval algorithm is optimal estimation (OE) (Rodgers, 2000). Therefore, DISAMAR provides full diagnostic information, including error covariance matrices and averaging kernels. It can deal with combined retrievals using information from different wavelength bands. OE is used for ozone profile retrieval by DISAMAR. The retrieved ozone profiles consist of the volume mixing ratio (VMR) of ozone at each altitude, obtained through spline interpolation or linear interpolation on the logarithm of the VMR given at a userdefined altitude grid. This makes conversion from VMR to Dobson units per layer and from VMR to number densities simple.
By using the newly developed DISMAS approach, which combines the principle of the DOAS retrieval method with OE, the number of wavelengths at which forward calculations have to be performed can be reduced significantly.
The forward model of DISAMAR provides simulated radiance (or sunnormalized radiance) and irradiance and the derivatives of the radiance with respect to the state vector elements. The semianalytical approach computes the forwardmode (tangentlinear) derivatives of the radiance spectrum. Here, attention is given to the calculation of these derivatives using a semianalytical approach because this important aspect for retrievals has not been published before.
3.1 Instrument model
The simulated radiance and irradiance can be convoluted separately with the instrument spectral response function (ISRF). The ISRF is defined in terms of a convolution operator. A Gaussian ISRF and a flattopped ISRF are included in DISAMAR. A tabulated ISRF, e.g. measured for a specific instrument, can also be used as an external file. During the convolution, the wavelength grid of the ISRF is interpolated to a highresolution wavelength grid (see Sect. 3.2.1 for the highresolution wavelength grid). The instrument features, such as signaltonoise ratio (SNR), stray light, and wavelength calibration errors, can be simulated using DISAMAR by specifying corresponding entries in a configuration file.
3.2 Atmospheric and surface properties
3.2.1 Gas absorptions
Gas absorptions include those from line absorbers and nonline absorbers. In the wavelength region considered here (270 to 2400 nm), line absorbers are H_{2}O, CO_{2}, CO, CH_{4}, and O_{2}. Nonline absorbers are O_{3}, SO_{2}, O_{2}O_{2}, BrO, HCHO, NO_{2}, and CHOCHO. The collision complexes, such as O_{2}O_{2} and N_{2}O_{2}, absorb light at UV, Vis, NIR, and SWIR (shortwave infrared) wavelengths. The amount of absorption by the collision complexes generally scales with the square of the pressure. In DISAMAR absorption by N_{2}O_{2} (in the spectral region of the O_{2} Aband) and the absorption by O_{2}O_{2} are simulated.
The standard database for lineabsorbing molecules is the HITRAN 2008 database (Rothman et al., 2009). Line parameters are read from the HITRAN database for a particular gas and a Voigt profile is used to calculate the absorption cross section. For H_{2}O, CO_{2}, CO, and CH_{4}, line mixing is ignored. For O_{2}, line mixing is taken into account using the model described in Tran and Hartmann (2008) for the O_{2} Aband (De Haan, 2012).
The gas absorptions from line absorbers have to be calculated linebyline at a highresolution wavelength grid, and then convolved with the ISRF. Gaussian division points are used to define the wavelength grid, which improves the accuracy of the integration during the convolution of the ISRF.
There are two steps to construct the highresolution wavelength grid: (1) determining wavelength intervals and (2) dividing each interval with proper Gaussian division points. The full width at half maximum (FWHM), the number of Gaussian division points for one FWHM (N_{0}), and the minimum number of Gaussian division points are specified in the configuration file.

We start from the shortest wavelength with an interval of one FWHM of the ISRF. If there is a strong absorption line in the interval, the boundary of the interval is set to the position of the strong line, making this interval smaller than one FWHM. The next interval starts from the position of the strong line with one FWHM interval again. This process is repeated until the end of the wavelength range.

If the interval is one FWHM, the interval is divided by the number of Gaussian points. If the interval is smaller than one FWHM, the number of Gaussian points is scaled with the size of the interval. For example, if the interval is half of one FWHM, the number of Gaussian points is 0.5×N_{0}. If the scaled number of Gaussian points is smaller than the minimum number of Gaussian points, the minimum number of Gaussian points is used. Note that the Gaussian quadrature weights and abscissae are determined for each interval.
Therefore, the wavelength grid is not equidistant: a finer wavelength grid is used for denser absorption lines, and a coarser grid is created if there are no absorption lines. Typically the number of Gaussian division points is between 3 and 30 per wavelength interval in the O_{2} Aband for a FWHM of 0.5 nm.
Figure 1 shows two 1 nm wide parts of the O_{2} Aband spectrum, from 757.8 to 758.8 and from 765.2 to 766.2 nm, simulated using DISAMAR. The used spectral resolution (wavelength difference between two grid points) is also shown in Fig. 1. In the spectral part from 757.8 to 758.8 nm, which is devoid of O_{2} absorption lines, the wavelength intervals are 0.4 nm, and every interval is divided into 40 Gaussian points (see Fig. 1a, b). In the spectral part from 765.2 to 766.2 nm, which has many O_{2} absorption lines, the intervals are smaller due to the absorption lines (see Fig. 1c, d).
Vertical profiles of absorbing gases are specified using the volume mixing ratio (VMR, in parts per million by volume, ppmv) at pressure grid levels. The volume mixing ratio has the advantage that it remains constant if the temperature changes. The number density of gas molecules is calculated internally using hydrostatic equilibrium. The pressure profile and temperature profile grid can be different for simulation and retrieval: for example, trace gas profile retrievals may be using 50 grid points for the simulation to have high forward model accuracy but only 12 or 18 grid points for the retrieval to represent the limited information content of the measurements.
3.2.2 Cloud and aerosol layers
Two different types of cloud and aerosol layers are distinguished in DISAMAR: an opaque Lambertian reflector and a layer of scattering particles. For the Lambertian type, the cloud and aerosol reflector is located at the top of a pressure interval. For the scattering particles layer, a Henyey–Greenstein phase function or a Mie phase function (or phase matrix in case of polarization) can be used for aerosols and clouds in DISAMAR. If the Henyey–Greenstein phase function is used, the wavelength dependence of the optical thickness is modelled using the Ångstrøm exponent and the absorption is controlled by the singlescattering albedo of the aerosol and cloud particles. A Mie phase function is provided in the form of expansion coefficients (De Rooij and Van der Stap, 1984) in an external file. It is assumed that in general a part of the pixel is cloudfree and that a part is covered by cloud. This part is controlled by the cloud fraction, which can have values in the interval [0, 1].
In DISAMAR the atmosphere is vertically divided into pressure intervals. The surface pressure (e.g. 1000 hPa) defines the lower boundary of the first interval. The intervals are further specified by the pressures of the top of the intervals. For example, if the pressure levels are specified at 700, 600, and 0.1 hPa, the atmosphere has three intervals: from the surface to 700 hPa, from 700 to 600 hPa, and from 600 to 0.1 hPa. The top of the atmosphere (TOA) is at 0.1 hPa. This division is suited to model a cloud layer with boundaries at 700 and 600 hPa. Internally, the pressure levels are translated into altitude levels assuming hydrostatic equilibrium. The formulae are provided in Appendix C. For this calculation, the temperatures at the pressure levels should be provided in the configuration file.
In order to perform radiative transfer calculations, each interval is divided into a number of homogeneous layers using Gaussian division points for the altitude. For example, the three intervals from the surface to TOA may be divided into layers using 12, 10, and 28 Gaussian division points, respectively. The standard atmospheres can be used as input for the pressure, temperature, and gas mixing ratio profiles.
Figure 2 illustrates the pressure intervals and the layers in the intervals used in DISAMAR. It is important to use a sufficient number of Gaussian division points for each pressure interval. Each layer is actually further divided into sublayers using Gaussian division points, with typically four sublayers per layer. In this way, an accurate integration over altitude can be performed in order to get the layer quantities. Averaging of the optical properties and radiation quantities at the top and bottom of a layer is not sufficiently accurate if the optical properties differ strongly between the top and bottom of the layer. The VMR profiles of gases are interpolated at the radiative transfer model grid. For optically thick clouds, we suggest using 15–50 Gaussian division points for the cloud layer, particularly if the cloud parameters for that layer are to be fitted. For optically thin clouds, 10 Gaussian division points should be sufficient for the cloud layer. We use Gaussian division points because it improves the accuracy in the integration over the altitude.
3.2.3 Surface albedo
The surface below the atmosphere is a Lambertian, i.e. isotropically reflecting, surface. In DISAMAR one can choose between a wavelengthindependent surface albedo or a surface albedo that has a polynomial wavelength dependence for each spectral band.
Suninduced fluorescence from vegetation due to photosynthesis can be observed in the nearinfrared spectrum (Joiner et al., 2011; Frankenberg et al., 2011). In order to be able to fit fluorescence in DISAMAR, a polynomial describing surface emission can be used in combination with a wavelengthdependent surface albedo.
3.3 Radiative transfer
3.3.1 Doubling–adding method
The doublingandadding method (also called adding method) is described in many books and papers on radiative transfer in planetary atmospheres (e.g. Van de Hulst, 1981; Hansen and Travis, 1974; Hovenier et al., 2004). Here we use the notation that follows De Haan et al. (1987) and Hovenier et al. (2004). We present the new formulae needed to calculate the derivatives; for the detailed formulae used in the doubling–adding method, we refer to De Haan et al. (1987) and Hovenier et al. (2004).
DISAMAR uses the doubling method for the calculation of the scattering properties of the individual atmospheric layers. The adding method is used to add different atmospheric layers, and thus it is used to construct the entire atmosphere. As discussed by De Haan et al. (1987) and Stammes et al. (1989), the doubling–adding method can be used to calculate the internal radiation field at the interface between layers in an efficient manner using principles of invariance. This approach can be extended to obtain the derivatives of the internal radiation field, in a linearization step, with respect to a parameter of the atmosphere–surface system. Because the adding formulae are needed in the derivation of the linearization, we provide them here as well.
Radiative transfer matrices are denoted in bold font, and they operate on incident radiation from right to left. A beam of polarized radiation is represented by the Stokes fourvector $\mathit{I}=[I,Q,U,V{]}^{T}$. Incident sunlight is the Stokes fourvector ${\mathit{E}}_{\mathrm{0}}=[{E}_{\mathrm{0}},\mathrm{0},\mathrm{0},\mathrm{0}{]}^{T}$, where E_{0} is the extraterrestrial solar irradiance. The elements of the 4×4 scattering matrix F are expanded into generalized spherical functions, where the expansion coefficients are calculated using the Mie computer program (De Rooij and Van der Stap, 1984). A generalized form of the addition theorem is used to calculate the Fourier coefficients of the 4×4 phase matrix, denoted by Z. Z is found by rotating F to the local meridian plane of incident and scattered light, respectively. Figure 3 illustrates the definition of the directions used in DISAMAR.
3.3.2 Radiation quantities for an optically thin layer
The doubling method is started with an optically thin layer, with an optical thickness of, e.g. $\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$. For a thin layer, the Fourier coefficients (indicated by m) of the reflection matrix R^{m} and the transmission matrix T^{m} of the layer illuminated at the top are given by
If the thin layer is illuminated at the bottom, the reflection and transmission matrix Fourier coefficients are given by
Here k_{sca} is the volume scattering coefficient of the layer, z is the altitude, and μ and μ^{′} are the cosines of the polar angles for scattered and incident light, respectively. The “*” indicates quantities with illumination at the bottom. ${\mathbf{Z}}^{m}(\pm \mathit{\mu},\pm {\mathit{\mu}}^{\prime})$ are the Fourier coefficients of the phase matrix of the layer. We assume that the scattering matrix does not depend on altitude within a homogeneous layer The unit of k_{sca} is per kilometre if z is in kilometres.
Attenuation of an incident beam of light due to extinction (sum of scattering and absorption) by the thin layer is given by the direct transmission matrix E_{thin}:
where the volume extinction coefficient is ${k}_{\mathrm{ext}}\left(z\right)={k}_{\mathrm{sca}}\left(z\right)+{k}_{\mathrm{abs}}\left(z\right)$, with k_{abs}(z) as the volume absorption coefficient. $\mathit{\delta}(\mathit{\mu}{\mathit{\mu}}^{\prime})$ is the Dirac delta function, and 1 is the unity matrix. Equation (5) follows from the Lambert–Beer attenuation law by keeping only the linear term of a Taylor series expansion.
3.3.3 Adding method formulae
We now consider adding two layers on top of each other. Let E_{k}, T_{k}, and R_{k} be the matrices representing the Fourier coefficients of the direct transmission, diffuse transmission, and reflection for illumination at the top of layer k, respectively. Here we omitted the Fourier index m for brevity. Similarly, let ${\mathbf{E}}_{k}^{*}$, ${\mathbf{T}}_{k}^{*}$, and ${\mathbf{R}}_{k}^{*}$ be these matrices for illumination at the bottom of layer k. The matrices for layer k+1 are denoted similarly. E_{k} is a diagonal matrix. For a planeparallel atmosphere, its elements are given by
where i is the index of the Gaussian μ point, $j=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}$ is the Stokes parameter index, and τ is the layer optical thickness.
The adding scheme for illumination at the top of the two layers is presented in Eqs. (7)–(12). The adding scheme for illumination at the bottom of the two layers is presented in Appendix B. The adding scheme is used for every Fourier coefficient; at the end of the process, the reflection and transmission matrices are obtained by summing up the Fourier coefficients.
The scheme of adding layer k+1 on top of layer k to create a combined layer is illustrated in Fig. 4. The adding method formulae are as follows.
Here ${\mathbf{E}}_{k,k+\mathrm{1}}$ is the attenuation of incident light by the combined layer. ${\mathbf{Q}}_{k,k+\mathrm{1}}$ is the sum of the repeated reflections between the layers for a unit amount of light incident on the interface between the layers. ${\mathbf{D}}_{k,k+\mathrm{1}}$ is the diffuse downward light at the interface between the layers, ${\mathbf{U}}_{k,k+\mathrm{1}}$ is the upward light at the interface between the layers, ${\mathbf{R}}_{k,k+\mathrm{1}}$ is the reflectance of the combined layer, and ${\mathbf{T}}_{k,k+\mathrm{1}}$ is the transmission of the combined layer.
The matrices R, T, E, Q, D, and U are socalled supermatrices. The elements of supermatrices contain Gaussian points for the polar angle with weights that are suited for numerical integration by multiplication of two supermatrices. For the definition of supermatrices, we refer to De Haan et al. (1987). The total number of polar angles, N_{μ}, is the number of Gaussian division points, N_{g}, plus one solar zenith angle and one viewing zenith angle, pertaining to the selected sunview geometry. The dimension of the supermatrix is (N_{μ}, N_{μ}) if polarization is ignored and (4N_{μ}, 4N_{μ}) if polarization is included using all four Stokes parameters. In the following a supermatrix is meant when we speak of a matrix.
The formula for repeated reflections at the interface, Eq. (8), can also be written as
where 1 is the unity matrix. Inserting Eq. (13) into Eq. (9) yields
Equations (13) and (14) are needed in the derivation of the derivatives, as will be discussed next.
3.3.4 Derivatives of the reflectance using the adding method
Retrieval algorithms that are designed to derive quantitative information about properties of the atmosphere from measured spectra require not only an atmospheric model but also derivatives of the reflectance with respect to the optical properties of the atmosphere, sometimes called Jacobians or weighting functions (see, e.g. Spurr, 2006; Hasekamp and Landgraf, 2005; Bai et al., 2020; Rodgers, 2000), to simulate reflectance spectra. Obviously, these derivatives can be calculated using numerical differentiation based on perturbation. However, numerical differentiation is a slow process, in particular when altitude profiles have to be determined and the atmosphere consists of 20–100 layers. Here we derive semianalytical expressions for the derivatives of the reflectance with respect to the optical properties of the atmospheric layers. The semianalytical derivatives are calculated at the interfaces of the layers, i.e. at levels. The derivatives for the layers themselves can then be calculated from the derivatives at the interfaces using interpolation and integration over the layers. Note that often the derivatives are calculated pertaining to optical properties of atmospheric layers, which differs from our approach.
The derivatives can be found by calculating the change in reflectance at topofatmosphere using the changes in atmospheric properties at each height, $\mathrm{d}\mathbf{R}/\mathrm{d}z$. In the adding scheme, such a change can be invoked by adding a thin layer at altitude z. As shown by Eqs. (1)–(5), the reflection and transmission of a thin layer can be written analytically as a function of the optical properties of the layer: k_{sca}, k_{abs}, and Z. Therefore, the computation of semianalytical derivatives should involve a thin layer. If one first calculates the reflection of two layers on top of each other, R_{top+bot}, and next calculates the reflection of the two layers with a thin layer in between, ${\mathbf{R}}_{\mathrm{top}+\mathrm{thin}+\mathrm{bot}}$, then $\mathrm{d}\mathbf{R}/\mathrm{d}z$ can be calculated from the following differential: (${\mathbf{R}}_{\mathrm{top}+\mathrm{thin}+\mathrm{bot}}{\mathbf{R}}_{\mathrm{top}+\mathrm{bot}})/\mathrm{d}z$.
We divide the atmosphere into two parts, a top part and a bottom part. Let R_{top}, ${\mathbf{R}}_{\mathrm{top}}^{*}$, T_{top}, ${\mathbf{T}}_{\mathrm{top}}^{*}$ be the Fourier coefficients of the reflection and transmission matrices of the top part; the quantities without an asterisk denote the properties for illumination at the top, and the quantities with an asterisk denote illumination at the bottom. Similarly, matrices with the subscript “bot” denote the properties of the lower part of the atmosphere.
Using the adding formulae, Eqs. (7)–(12) and keeping only constant terms and terms linear in dz, one can derive the derivative $\mathrm{d}\mathbf{R}/\mathrm{d}z$ in a concise form as follows:
where Q and Q^{*} are given by
and W(z) is given by
Here $\mathbf{1}\mathbf{/}\mathit{\mu}$ is the supermatrix with the element $\mathrm{1}/\mathit{\mu}$ on the diagonal, representing the slant path for extinction. $\mathbf{Z}{{}^{\prime}}_{\pm \pm}$ are abbreviated notations of $\mathbf{Z}{}^{\prime}(\pm \mathit{\mu},\pm {\mathit{\mu}}^{\prime})$, the supermatrices of the Fourier coefficients of the phase matrix. The full derivation of Eq. (15) is given in Appendix A.
Inserting Eq. (18) into Eq. (15) yields the separate terms of the derivative:
Since we want to calculate the derivative $\mathrm{d}\mathbf{R}/\mathrm{d}z$ with respect to the scattering and absorption properties of the atmosphere at altitude z, Eq. (19) has to be expressed in the internal fields D and U at altitude z. Therefore, in Eqs. (20)–(23) the internal fields are expressed using (1+Q) or ($\mathbf{1}+{\mathbf{Q}}^{*}$) and R, T, and E. The relations are as follows.
Here $\stackrel{\mathrm{\u0303}}{\phantom{\rule{0.25em}{0ex}}}$ indicates the transpose of a supermatrix, and ${\mathbf{\Delta}}_{\mathrm{3}}=\mathrm{diag}(\mathrm{1},\mathrm{1},\mathrm{1},\mathrm{1})$ is a diagonal matrix that occurs due to polarization. By replacing the terms related to (1+Q) and ($\mathbf{1}+{\mathbf{Q}}^{*}$) in Eq. (19) with the relations in Eqs. (20)–(23), we obtain the derivative in terms of the internal radiation fields:
First, replace k_{ext} with k_{sca}+k_{abs} in Eq. (24) and separate the terms containing k_{abs} and k_{sca}. Following this, calculate the partial derivative of $\mathrm{d}\mathbf{R}/\mathrm{d}z$ with respect to k_{abs} and k_{sca} using Eq. (24). This gives the partial derivative of the reflectance with respect to the volume absorption coefficient (Eq. 25) and the partial derivative of the reflectance with respect to the volume scattering coefficient (Eq. 26). These equations are as follows:
The above relation, Eq. (26), is a generic expression for the derivative with respect to the volume scattering coefficient. However, since it contains the phase matrix, it makes a difference whether the change in the volume scattering coefficient is due to a change in the number of scattering molecules in the thin layer or due to a change in the number of aerosol particles. If the number of scattering molecules is changing, the phase matrix of the molecules has to be used in Eq. (26), whereas the phase matrix of the aerosol particles has to be used if the number of aerosol particles is changing. In general one has to distinguish between molecules, aerosols, and cloud particles.
3.3.5 Derivatives of the reflectance with respect to specific parameters
If polarization is ignored, the partial derivative of the reflectance with respect to some specific quantity x(z) (e.g. ozone number density or aerosol number density) can be expressed as follows:
where R is the (1, 1) element of R and ${\mathit{\alpha}}_{\mathrm{1}}^{l}$ is the phase matrix expansion coefficient for unpolarized light. In DISAMAR the third term on the righthand side is ignored, so it is assumed that the phase matrix of the particles does not change as function of z. The derivatives for specific parameters can be calculated using these basic derivatives (De Haan, 2012). Here we give two examples of the use of the semianalytical derivatives presented in the previous section.
The first example is the altituderesolved air mass factor (also called block air mass factor) at wavelength λ, m(z,λ). The altituderesolved air mass factor is used in DOAS retrievals and represents the relative reduction in the reflectance when a unit amount of absorption is added to the atmosphere in a thin layer located between z and z+dz. So it is the partial derivative of the reflectance due to the absorber at z divided by the reflectance. Writing the reflectance of the atmosphere without the extra absorption as R(λ,0), the altituderesolved air mass factor is defined as follows:
The partial derivative on the righthand side can now be calculated from Eq. (25).
The second example is the derivative of the reflectance to the altitude of a Lambertian cloud below an absorbing and scattering gas layer:
Here the cloud is approximated as a Lambertian surface, without any transmission. When the altitude of the Lambertian cloud increases, the amount of gas above the cloud is reduced, and thus we get a minus sign in Eq. (29). The amount of scattering gas that is removed is proportional to the volume scattering coefficient (k_{gas,sca}), whereas the amount of absorbing gas that is removed is proportional to the volume absorption coefficient (k_{gas,abs}). The partial derivatives on the righthand side of Eq. (29) can now be calculated from Eqs. (25) and (26).
As noted before, the derivatives can also be calculated numerically. Thus, $\mathrm{d}R(\mathit{\lambda},\mathit{x})/\mathrm{d}{x}_{k}$, with x_{k} element k of the state vector x, is written as the following numerical differential:
Here Δx=0, except for the element k, where Δx=Δx_{k}. By comparing the analytical derivatives in DISAMAR with the numerical derivatives, we can validate the former.
3.3.6 Layerbased orders of scattering method (LABOS)
In the adding method as implemented in DISAMAR, the reflection matrix is calculated for a set of solar directions equal to the number of Gaussian division points used for integration over the polar angle plus at least two additional directions, namely the actual viewing direction and the actual solar position. Thus, for satellite retrievals the adding method provides more results than required.
The layerbased orders of scattering method (LABOS) has been developed to obtain the reflection matrix for the actual solar position or, when derivatives are required, for the actual solar position and a solar position corresponding to the viewing direction. The reflection and transmission properties of the individual homogeneous layers are still calculated with the doubling method; only the adding of different layers and the subsequent calculation of the internal field is replaced by the successive orders of scattering method (see Eqs. 35–38). Here an order of scattering represents scattering by an atmospheric layer. This differs from the classical method of successive orders of scattering where the scattering element is a volume element of the atmosphere instead of a layer (Lenoble et al., 2007; Min and Duan, 2004). In the adding method, one deals with matrix–matrix multiplications, whereas in LABOS one deals with matrix–vector multiplications. However, in LABOS the calculations have to be repeated for the different orders of scattering. Hence, for an optically thin atmosphere with only a few orders of scattering, LABOS requires less calculation time, while adding is more efficient for optically thick atmospheres, e.g. those containing clouds.
In LABOS we deal with the upward and downward internal radiation fields, i.e. U and D. In order to determine the derivatives using reciprocity, we need to consider two directions for the incident light, one corresponding to the actual solar position and one corresponding to the viewing direction. Therefore, U and D are (4N, 2) matrices, with the first column corresponding to incident light for the viewing direction, and the second column corresponding to the solar direction, where the number of polar angles ${N}_{\mathit{\mu}}={N}_{\mathrm{g}}+\mathrm{2}$. If the incident light is polarized, we need (4N_{μ}, 8) matrices. In order to calculate the total internal fields U and D (see Eqs. 35–38), we first need to calculate local internal fields U^{local} and D^{local} (see Eqs. 31–32).
Let R_{k} and T_{k} be the supermatrices for the reflection and transmission of layer k ($k=\mathrm{0},\mathrm{1},\mathrm{\dots}K$). If k=0, the reflection matrix is that of the surface, and the transmission matrix vanishes. Further, let μ and μ_{0} be the cosines of the viewing and solar zenith angle, respectively. In order to evaluate different orders of scattering, we first consider the socalled local radiation field, which is due to the reflection by the layer just below the level considered and the transmission by the layer just above this level. Figure 5 illustrates the local upward and downward radiation fields at the interface k for the firstorder scattering (n=1).
The local upward firstorder radiation at level k is given by
where τ_{k} is the vertical optical depth of level k, measured from the top of the atmosphere. The two columns of matrix ${\mathbf{U}}_{k}^{\mathrm{local},\mathrm{1}}$ correspond to the two columns of supermatrix R_{k} having indices 4N_{g}+1 and $\mathrm{4}({N}_{\mathrm{g}}+\mathrm{1})+\mathrm{1}$, respectively, if the dimension of the Stokes vector is 4. The firstorder local downward radiation field at level k is
Here the two columns of ${\mathbf{D}}_{k}^{\mathrm{local},\mathrm{1}}$ correspond to the two columns of supermatrix T_{k+1} with indices 4N_{g}+1 and 4N_{g}+5. The superscript 1 denotes the firstorder scattering for the layers.
From the local radiation fields, we can now derive the total radiation fields. The firstorderscattered upward radiation at level k, ${\mathbf{U}}_{k}^{\mathrm{1}}$, is calculated as the sum of the contributions of the local upward radiation from the ground up to level k. The downward radiation at level k, ${\mathbf{D}}_{k}^{\mathrm{1}}$, is calculated as the sum of the contributions of the local downward radiation from level k up to level K (TOA):
where the attenuation accounts for the direct transmission of light from level l to level k. We note that Eqs. (33) and (34) can be rewritten in terms of the following recurrence relations, which are more efficient to evaluate.
The same recurrence relations hold for the higher orders of scattering because the direct transmission of radiation to other levels is the same for each order of scattering.
The following equations are used to determine the local radiation fields for the higher orders of scattering, where n indicates the order of scattering:
where ${\mathbf{R}}_{k+\mathrm{1}}^{*}$ is the reflection of layer k+1 for illumination at its bottom. Since the layers are homogeneous, we can use symmetry relations to calculate ${\mathbf{R}}_{k+\mathrm{1}}^{*}$ and ${\mathbf{T}}_{k}^{*}$ from R_{k+1} and T_{k}, respectively (Hovenier et al., 2004). After the local radiation fields are calculated using Eqs. (39) and (40), the recurrence relations are applied for the transfer to other levels to provide the radiation fields ${\mathbf{U}}_{k}^{n+\mathrm{1}}$ and ${\mathbf{D}}_{k}^{n+\mathrm{1}}$ for scattering order n+1.
One can continue calculating the orders of scattering and summing them up until the contribution of higher orders can be ignored. At the end of the calculation, the internal radiation fields at all levels are known. The Fourier coefficients of the upward radiation at the top of the atmosphere provide the reflection, and the internal radiation fields are used to calculate the derivatives. In DISAMAR users can choose to use either the doubling–adding method or the doubling–LABOS method.
3.3.7 Reflection calculated by integration over the source function
The upward radiation at the top of the atmosphere is calculated by assuming that the atmosphere is consisting of homogeneous layers, both in the doubling–adding method and in the doubling–LABOS method. The calculation will converge to the correct value if many homogeneous layers are used to represent the true atmospheric scattering and absorption profile. In DISAMAR the reflection can also be calculated from numerical integration over the source function J to mitigate errors due to layering that is too coarse. The reason is that by using the source function the assumption of homogeneous layers is not used for single scattering but only for multiple scattering. For each azimuthal Fourier term the reflection is then given by
where the surface is at z=0. Since the source function for the surface is a Dirac delta function in z, its contribution to the reflection occurs as a separate term in Eq. (41). The source function for the atmosphere is
$\mathbf{U}(z,{\mathit{\mu}}^{\prime},{\mathit{\mu}}_{\mathrm{0}})$ and $\mathbf{D}(z,{\mathit{\mu}}^{\prime},{\mathit{\mu}}_{\mathrm{0}})$ are the internal fields at altitude z that are known after doubling–adding or doubling–LABOS calculations. Hence, J can be easily calculated at those altitudes where the internal radiation fields are known.
Three retrieval methods are implemented in DISAMAR: the optimal estimation (OE), differential optical absorption spectroscopy (DOAS), and differential and smooth absorption separated (DISMAS) methods.
The OE method is implemented in DISAMAR following Rodgers (2000). The OE method can be used for various retrievals implemented in DISAMAR, especially ozone profile retrieval (Kroon et al., 2011) and aerosol layer height retrieval (Sanders et al., 2015; Nanda et al., 2018). The OE output includes the retrieved parameters and their error matrices, correlations, gain matrices, and averaging kernels.
The DOAS method is described in the literature (Platt and Stutz, 2008). Below we describe the DISMAS method, which is related to the DOAS method, in detail.
4.1 DISMAS method
The DISMAS method uses the principle of DOAS for the forward simulation of spectra and the OE method for the retrieval. Therefore, the DISMAS method reduces the forward simulation time, making it faster than a full forward simulation, but still provides the full error estimation as in OE. DISMAS can be used for the retrieval of total columns of weakly absorbing gases, like O_{3} (in the visible wavelength range), NO_{2}, SO_{2}, BrO, and CH_{2}O, but not for strong absorbers, such as O_{2}, H_{2}O, CH_{4}, CO, and O_{3} (in the UV). In DISMAS, gas absorption spectra are separated into a smooth part and a differential part, which is similar to DOAS. As an illustration, Fig. 6 shows an NO_{2} absorption spectrum with a smooth and a differential absorption part.
For the reflectance spectrum around a weak atmospheric absorption line, e.g. from NO_{2}, one can use the Lambert–Beer attenuation law as follows:
where R(λ,0) is the reflectance of the atmosphere without the absorbing gas, M(λ) is the air mass factor (AMF) of the total atmosphere (unitless), ${\mathit{\sigma}}_{\mathrm{a}}^{\mathrm{abs}}\left(\mathit{\lambda}\right)$ is the altitudeaveraged (subscript a) absorption cross section (in ${\mathrm{cm}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{molec}.}^{\mathrm{1}}$), and N is the vertical column density of the absorbing gas (in $\mathrm{molec}.\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{\mathrm{2}}$). Introducing the smooth absorption cross section, ${\mathit{\sigma}}_{\mathrm{a}}^{\mathrm{sm}}$, and the differential absorption cross section, ${\mathit{\sigma}}_{\mathrm{a}}^{\mathrm{dif}}$, one obtains
Here $R(\mathit{\lambda},{\mathit{\sigma}}_{\mathrm{a}}^{\mathrm{sm}}(\mathit{\lambda}\left)\right)$ and M(λ) are smooth functions of wavelength, which are calculated with a radiative transfer model at a few wavelengths and are fitted as a function of wavelength to be interpolated to a highresolution wavelength grid. The differential term
varies strongly with wavelength and has to be evaluated at a highresolution wavelength grid. However, this analytical formula can be calculated very quickly without the need to solve the radiative transfer model. The smooth functions $R(\mathit{\lambda},{\mathit{\sigma}}_{\mathrm{a}}^{\mathrm{sm}}(\mathit{\lambda}\left)\right)$ and M(λ) are calculated at those wavelengths, λ_{p}, where the differential absorption cross sections are zero because at those points (see Eq. 44)
Hence, radiative transfer calculations for a small number of wavelengths λ_{p} are sufficient to determine the reflectance spectrum at a high spectral resolution grid.
To simulate a measured spectrum in DISMAS, the modelled reflectance spectrum on the highresolution grid is multiplied with a highresolution solar irradiance spectrum and convoluted with the instrument's spectral response function (or slit function). This provides a simulated radiance spectrum of the instrument that can be compared to an actually measured radiance spectrum. In practice, the simulated radiance spectrum is divided by the simulated irradiance spectrum to eliminate common errors, and the logarithm of the reflectance is used for fitting.
The air mass factor M(λ) can be calculated from m(z,λ), which is known from the radiative transfer calculations (see Eq. 28):
where n(z) is the number density profile of the absorbing gas. To start the optimal estimation process in DISMAS, m(z,λ) is initialized using a geometrical air mass factor. In subsequent iteration steps, the value of m(z,λ) is updated by radiative transfer calculations for the new atmospheric composition state vector.
If there is more than one absorbing trace gas in a fit window, one can use the slant absorption optical thickness, ${\mathit{\tau}}_{\mathrm{slant}}^{\mathrm{abs}}\left(\mathit{\lambda}\right)$, instead of the absorption cross section to separate the smooth and differential parts of the absorption spectrum.
For fitting the modelled reflectance to the measured reflectance and determining the errors in the retrieved parameter values, the derivatives are needed. Here we give the derivatives with respect to the total column of a trace gas and to the cloud fraction for a partly cloudy pixel.
For a partly cloudy pixel with K different trace gases, the reflectance is written as follows:
where c is the cloud fraction, subscript “sm” refers to the smooth reflectance, and superscripts “clr” and “cld” denote the clear and cloudy part of the pixel, respectively. The derivative with respect to the cloud fraction is
The derivative with respect to the total column of trace gas k, N_{k}, becomes
where we assume that the air mass factor does not change significantly when the column changes, which is accurate for weak absorption.
5.1 Comparison with the Doubling–Adding KNMI (DAK) model
As an example, we show a comparison of DISAMAR with the DAK radiative transfer model (De Haan et al., 1987; Stammes et al., 1989) for the simulated reflectance and linear polarization spectrum of a Rayleigh scattering atmosphere in the ozone Huggins band region between 325 and 335 nm; see Fig. 7. In DISAMAR we used the doublingadding method in order to get the best agreement with DAK. The atmospheric temperature, pressure, and ozone mixing ratio profiles are taken from the midlatitude summer atmosphere (Anderson et al., 1986). Since DISAMAR uses the hydrostatic relation to calculate altitude with the input pressure and temperature profiles, the altitude in DAK is corrected for the hydrostatic relation. This correction is important to obtain the same Rayleigh scattering optical thickness in this wavelength range. The total ozone column is 335.66 DU in DAK and DISAMAR. The ozone cross section was taken from Bass and Paur (1985) because these data are available in both DAK and DISAMAR.
Figure 7 shows the simulated reflectance, degree of linear polarization, and the corresponding relative differences, for the viewing and solar geometry μ=1, μ_{0}=0.5. DISAMAR and DAK agree very well for the simulated reflectance and degree of linear polarization, with absolute differences of less than 0.01 %. DISAMAR and DAK use different internal wavelength grids, altitude grids, and interpolation methods, which may explain the remaining small differences.
5.2 Derivatives
The accuracy of the semianalytical formulae for the derivatives based on the adding method is checked here by comparison with the numerically calculated derivatives. We show an example of the derivative of the cloud height, $\mathrm{d}R/\mathrm{d}z$, in the O_{2} Aband for a Lambertian cloud. In the calculation, the atmospheric profile is midlatitude summer, the cloud fraction is set to 1, and the cloud albedo set to 0.8. The viewing and solar geometry is μ=1, μ_{0}=0.8660. The spectrum of the derivative is calculated at a highresolution wavelength grid and convolved with a Gaussian ISRF with a FWHM of 0.4 nm. As shown in Fig. 8, the semianalytical and numerical derivatives are very similar: the mean relative difference is −0.13 %. At 758–759 nm, the O_{2} absorption is very small, about 100 times smaller than in the absorption band close to 760–761 nm. The spectrum of $\mathrm{d}R/\mathrm{d}z$ clearly shows that the cloud height information comes from the absorption band. If the cloud fraction is smaller than 1, the derivative is also smaller because the radiance is smaller. Although the relative difference between the semianalytical and numerical derivatives behaves steeply in the continuum at 758–759 nm, it is still within 0.2 % over the entire absorption band. Therefore, the semianalytical derivatives are considered accurate for cloud height retrievals.
The altituderesolved AMF, m(λ,z), is used in the DOAS algorithm to convert the slant column density of a gas to the vertical column density (see Sect. 3.3.5). The altituderesolved AMF was calculated using Eq. (28), and the semianalytical partial derivatives are taken from DISAMAR. It can also be calculated numerically, but that method is slower and requires a finer altitude grid to get accurate results. Figure 9 shows the altituderesolved AMF of NO_{2} calculated using the semianalytical derivatives for a cloudy scene and a clearsky scene. Here the atmospheric profile and NO_{2} mixing ratio are taken from the midlatitude summer atmospheric profile. A scattering cloud layer is specified between 800 and 700 hPa with a Henyey–Greenstein phase function, asymmetry parameter 0.85, and single scattering albedo 1.0. The cloud optical thickness is 10. The surface albedo is 0.05. The viewing and solar geometry is μ=1, μ_{0}=0.8660.
Figure 9 shows that inside the cloud the AMF has a peak close to the cloud top due to the enhancement from multiple scattering. Above the cloud top, the AMF is larger than the clearsky AMF due to the bright cloud layer. Below the cloud top, the AMF decreases and becomes smaller than the clearsky AMF. At TOA, the AMF for both scenes approach the geometric AMF, $\mathrm{1}/\mathit{\mu}+\mathrm{1}/{\mathit{\mu}}_{\mathrm{0}}$, which is 2.155 for this case. These are typical characteristics of the altituderesolved AMF of NO_{2}.
5.3 Comparison of DISMAS and full OE for NO_{2} retrievals
The DISAMAR code can simulate the retrieval of total columns of trace gases, ozone profile, cloud layer height, aerosol layer height, cloud (or aerosol) optical thickness, surface albedo, and surface pressure. Some OE retrieval applications were mentioned in the introduction of Sect. 4. Because the DISMAS retrieval approach is new and not used operationally, we show the DISMAS results for NO_{2} column retrieval in a cloudy scene here as an example. The scene in the simulation is constructed as follows: the atmosphere is a typical European polluted model, the cloud is a Lambertian cloud with a fixed albedo of 0.8 located at 4 km altitude, and the geometry is nadir view with solar zenith angle 60^{∘}. In the retrieval, the fitting window is 420–450 nm. Radiative transfer calculations are performed at four wavelengths for DISMAS and at about 900 wavelengths for OE. Thus, DISMAS enables a large reduction of the number of radiative transfer calculations needed for spectral fitting. Table 1 shows the comparison between DISMAS and full OE for the NO_{2} retrievals for two cases. In Case 1, three parameters are retrieved: NO_{2} total column density, surface albedo, and cloud fraction. In Case 2, the cloud fraction is fixed, and only the NO_{2} total column density and the surface albedo are retrieved.
The parameters were accurately retrieved by both methods. The error estimates for the retrieved surface albedo and cloud fraction differ only little, but DISMAS underestimates the error in the retrieved NO_{2} column.
If the cloud fraction is fixed to its true value (0.20) so that the NO_{2} column and the surface albedo are the only fit parameters, DISMAS and OE fully agree for the errors in the NO_{2} column, while the surface albedo values are also the same.
From our results, we find that DISMAS and OE are functionally equivalent for weak absorption, except that the error in the retrieved column may differ in some cases. The main advantage is that DISMAS is much faster than full OE, which might make operational use of DISMAS attractive. This would make it possible to derive the surface albedo and the cloud fraction directly from the measured radiance in the NO_{2} window, thereby reducing the bias of the retrieved NO_{2} column significantly.
The aim of the DISAMAR code is to determine instrument specifications and analyse methods for retrievals of atmospheric composition from satellite observations in the range of 270–2400 nm. We have described the principles of the DISAMAR radiative transfer model and its innovative retrieval capabilities. The novel features in DISAMAR are the semianalytical derivatives of the reflectance using the internal radiation field, the layerbased orders of scattering method LABOS, which speeds up the calculations, and the DISMAS retrieval method, which reduces the number of radiative transfer calculations in DOAStype retrievals. The semianalytical derivatives are based on the linearization of the adding method, which is detailed in Appendix A. Other features of DISAMAR that improve the accuracy of forward modelling are the Gaussian grids for height and wavelength and the integration over the source function to obtain the reflectance.
DISAMAR can be used as an accurate radiative transfer model for simulations and as a retrieval algorithm for several applications. We have given a few examples of applications, namely forward modelling of the polarized reflectance in the UV, the derivative of cloud height in the O_{2} Aband, and an application of DISMAS for NO_{2} retrieval in a cloudy scene. It is shown that for retrieval of total columns of weak absorbers, DISMAS has a similar functionality to optimal estimation but is much faster.
In the past few years DISAMAR has been used extensively in trace gas, aerosol, and cloud retrieval studies for OMI, TROPOMI, Sentinel4, and Sentinel5. DISAMAR is currently being used for operational retrievals of ozone profiles from UV spectra for OMI (Kroon et al., 2011) and TROPOMI (Veefkind et al., 2021) and aerosol layer height retrieval from the O_{2} Aband for TROPOMI (Nanda et al., 2019). More applications of DISAMAR are described in the DISAMAR algorithm document (De Haan, 2012), and can be found in the peerreviewed literature; for example, on the role of aerosols in retrievals of NO_{2}; see Castellanos et al. (2015) and Chimot et al. (2016, 2017).
Currently, DISAMAR is still under development and has certain limitations: Raman scattering in water (Vasilkov et al., 2002) is ignored, remote sensing from the ground and from aircraft is not supported, and limb observations are not supported. Furthermore, the radiative transfer part assumes a planeparallel atmosphere, albeit with a correction for the curvature of the atmosphere for incident sunlight (pseudospherical correction; see Appendix B). This may become inaccurate for observations far from the nadir (viewing zenith angle $>{\mathrm{60}}^{\circ}$), in particular for stratospheric gases (Spurr, 2002). In the coming years DISAMAR will be further extended and used to develop optimized retrieval algorithms, including those using neural networks (Nanda et al., 2019).
Here we derive the formulae for the derivative of the reflection, $\mathrm{d}\mathbf{R}/\mathrm{d}z$, given by Eq. (15) in Sect. 3.3.4, by means of linearization of the adding equations. To that end, we divide the atmosphere into two parts, a top part and a bottom part. Let R_{top}, ${\mathbf{R}}_{\mathrm{top}}^{*}$, T_{top}, and ${\mathbf{T}}_{\mathrm{top}}^{*}$ be the Fourier coefficients of the reflection and transmission matrices of the top part. The quantities without an asterisk denote the properties for illumination at its top, while the asterisk denotes illumination at its bottom. Analogously, matrices with the subscript “bot” denote the properties of the lower partial atmosphere. We start with the adding equations Eqs. (11), (10), and (9) and change subscript k+1 to “top” and subscript k to “bot”:
Inserting Eqs. (A2) and (A3) into Eq. (A1), we obtain the adding equation for the reflection at the top in a concise form:
where
We now start the linearization process by adding a thin layer on top of the lower partial atmosphere. Using Eq. (A4) and writing “thin” instead of “top” yields:
The radiation fields for a thin layer have already been given by Eqs. (1)–(5) in Sect. 3.3.2. We observe that R_{thin}, T_{thin}, ${\mathbf{R}}_{\mathrm{thin}}^{*}$, ${\mathbf{T}}_{\mathrm{thin}}^{*}$, and (1−E_{thin}) are proportional to dz. In the linearization process, we keep only terms that are constant or linear in dz. We note that only the first term of Eq. (A7) has to be kept because the second term is ∝(dz)^{2}. Inserting Eq. (A7) into Eq. (A6) and keeping only terms that are constant or linear in dz, the derivation steps towards the linearized form of Eq. (A6), given by Eq. (A8), are
We now add the upper partial atmosphere to this combined (“thin+bot”) lower layer and derive the linearized expression for ${\mathbf{R}}_{\mathrm{top}+\mathrm{thin}+\mathrm{bot}}$ given by Eq. (A9.). We start again from Eq. (A4):
Subtracting these two equations we get the difference in R due to the thin layer:
Let us now focus on the term between the square brackets in Eq. (A11). The linearization steps, in which only constant terms and terms linear in dz are kept, are
In step 3, ${\mathbf{Q}}_{\mathrm{top}+\mathrm{thin}+\mathrm{bot}}$ and Q_{top+bot} are replaced by the first term of their expansions; in step 7, the last term is removed because it depends on (dz)^{2}; in step 8, ${\mathbf{R}}_{\mathrm{bot}}{\mathbf{R}}_{\mathrm{top}}^{*}$ is approximated by ${\mathbf{Q}}_{\mathrm{top}+\mathrm{bot}}^{*}$. Because ${\mathbf{R}}_{\mathrm{thin}+\mathrm{bot}}{\mathbf{R}}_{\mathrm{bot}}$ is small, step 8 only introduces a small error. In order to get a compact form, in step 9 a small term, ${\mathbf{Q}}_{\mathrm{top}+\mathrm{bot}}^{*}({\mathbf{R}}_{\mathrm{thin}+\mathrm{bot}}{\mathbf{R}}_{\mathrm{bot}}){\mathbf{Q}}_{\mathrm{top}+\mathrm{bot}}$, is added.
Inserting Eq. (A12) into Eq. (A11), we obtain the linearized form of the difference dR:
Here the term $({\mathbf{R}}_{\mathrm{thin}+\mathrm{bot}}{\mathbf{R}}_{\mathrm{bot}})$ contains the optical properties of the thin layer: k_{sca}, k_{ext} and Z. This term can be written in terms of these optical properties by using Eq. (A8) and inserting the thinlayer formulae Eqs. (1)–(5).
For the definition of the phase matrix Z^{′} with special supermatrix weights we refer to De Haan, 2012; De Haan et al., 1987. If we now insert this relation, Eq. (A14), into Eq. (A13), we get the following end result:
This formula, Eq. (A15), is the derivative equation given in Eq. (15).
The pseudospherical correction in DISAMAR corrects for the curvature of the atmosphere for incident sunlight. For the other light paths the atmosphere remains planeparallel. In order to calculate the pseudospherical correction, the adding scheme for illumination of the layers from below is implemented. These equations are as follows:
We assume that the individual layers remain planeparallel; the only correction that is needed for the reflectance is to change the attenuation terms for incident sunlight. For the adding method it means replacing E_{k+1}, when it occurs on the righthand side of a term, with a new attenuation term. The pseudospherical correction method follows Spurr (2002), and thus the details are not given here.
The partial derivative of the reflectance to the absorption coefficient, Eq. (25) from Sect. 3.3.4, can be rewritten as
With the pseudospherical correction, the last term at the righthand side of Eq. (B6) should be replaced by Eq. (B9). Equations (B7)–(B9) present the derivation of the pseudospherical correction:
where R_{E} is the radius of the Earth in kilometres and z_{k} and z_{l} are the altitudes (in kilometres) of levels k and l, respectively. Here we use the local value of U to distinguish between the different levels below level k. The factor $\mathrm{exp}\left(\right({\mathit{\tau}}_{k}{\mathit{\tau}}_{l})/{\mathit{\mu}}_{\mathrm{0}})$ occurs in Eq. (B7) because we used the transpose of U; this is also the reason for the order of the arguments for ${\mathbf{U}}_{l}^{\mathrm{local},n}$. The sum over n is the sum over the orders of scattering. For a planeparallel atmosphere, Eq. (B7) can be written as Eq. (B8). For a spherical atmosphere, $\mathrm{exp}({\mathit{\tau}}_{l}/{\mathit{\mu}}_{\mathrm{0}})$ should be replaced by $\mathrm{exp}({\mathit{\tau}}_{l}^{\mathrm{slant}}({\mathit{\mu}}_{\mathrm{0}}\left)\right)$, noted as ${\mathbf{E}}_{\mathrm{top}}^{\mathrm{sph}}$ in Eq. (B10), and the factor $\mathrm{1}/{\mathit{\mu}}_{\mathrm{0}}$ should be replaced by $\mathrm{1}/\mathrm{cos}{\mathit{\theta}}_{\mathrm{0}}^{\prime}$, where $\mathrm{1}/\mathrm{cos}{\mathit{\theta}}_{\mathrm{0}}^{\prime}=\frac{({R}_{\mathrm{E}}+{z}_{k})}{\sqrt{({R}_{\mathrm{E}}+{z}_{k}{)}^{\mathrm{2}}({R}_{\mathrm{E}}+{z}_{l}{)}^{\mathrm{2}}(\mathrm{1}{\mathit{\mu}}_{\mathrm{0}}^{\mathrm{2}})}}$.
The partial derivative of the reflectance to the scattering coefficient, Eq. (26) from Sect. 3.3.4, with the pseudospherical correction becomes
This completes the formulae for the pseudospherical correction of the derivatives.
The pressure grid in DISAMAR is converted to an altitude grid using the hydrostatic equation (assuming hydrostatic equilibrium):
The scale height H(z,l) [m] is given by
where R=8.31451 is the universal gas constant [J K^{−1} mol^{−1}], T(z) is the temperature [K], g(z,l) is the gravitational acceleration [m s^{−2}], $M=\mathrm{28.964}\times {\mathrm{10}}^{\mathrm{3}}$ is the mean molecular mass of dry air [kg mol^{−1}], and l is the geographical latitude. The altitude at the surface is set to the altitude with respect to the geoid (mean sea level of the oceans).
An approximate expression for the gravitational acceleration as a function of the altitude above the mean sea level, z [m], and the geodetic latitude, l_{g}, is given by Cox (2000):
with
The integration over altitude is given by
where p_{1}=p(z_{1}) and p_{2}=p(z_{2}). In order to calculate the altitude grid, we start at the surface with altitude z_{s} and pressure p_{s}. Using f(z)=1, we get
Assuming $g(z,l)=g({z}_{s},l)$, using Eq. (C2), one first calculates H(ln (p),l), then evaluates Eq. (C6) for the pressure grid to get an estimation of the altitude grid. This procedure is repeated using the correct values of g(z,l) (calculated using Eqs. C3–C4) to get a new estimation of the altitude grid. If the new altitude grid is close to the previous altitude grid, the iteration is stopped. One or two iterations are needed to get the final altitude grid. Please note that the altitude grid is updated whenever the temperature profile changes.
The DISAMAR code v4.1.5 and data are available from the following Zenodo DOI: https://doi.org/10.5281/zenodo.6304984 (De Haan, 2022). The DISAMAR code is also available on GitLab: https://gitlab.com/KNMIOSS/disamar (last access: 2 September 2022).
JFdH developed the DISAMAR codes and wrote the algorithm document and manual. PW wrote the paper based on the algorithm document, performed the calculations, and made the figures. MS and JPV contributed to the DISAMAR code. MS and PW maintain the code. PS contributed to the writing of the paper. All coauthors commented on the paper.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Over the years many people have contributed to the testing of the code and development of the applications. We want to thank Mark ter Linden (Science and Technology B.V., Delft, the Netherlands) for building the software interface around DISAMAR.
The development of DISAMAR was funded by the Netherlands Space Office (NSO), under the OMI and TROPOMI science contracts, and by ESA under various future mission contracts.
This paper was edited by Volker Grewe and reviewed by Minzheng Duan, Patrick Stegmann, and Pengwang Zhai.
Anderson, G. P., Clough, S. A., Kneizys, F. X., Chetwynd, J. H., and Shettle, E. P.: AFGL (Air Force Geophysical Laboratory) atmospheric constituent profiles (0120 km). Environmental research papers, 19860515, Air Force Geophysics Lab., Hanscom AFB, MA (USA), https://apps.dtic.mil/sti/pdfs/ADA175173.pdf (last access: 2 September 2022), 1986. a
Bai, W., Zhang, P., Zhang, W., Li, J., Ma, G., Qi, C., and Liu, H.: Jacobian matrix for nearinfrared remote sensing based on vector radiative transfer model, Science China Earth Sciences, 63, 1353–1365, https://doi.org/10.1007/s1143001995867, 2020. a
Bass, A. M. and Paur, R. J.: The ultraviolet crosssections of ozone: I. The measurements II. Results and temperature dependence, in atmospheric ozone, in: Proceedings of the Quadrennial Ozone Symposium, edited by: Zerefos, C. and Ghazi, A., Halkidiki, Greece, 3–7 September 1984, Dordrecht, Reidel, 606–616, https://doi.org/10.1007/9789400953130_120, 1985. a
Castellanos, P., Boersma, K. F., Torres, O., and de Haan, J. F.: OMI tropospheric NO2 air mass factors over South America: effects of biomass burning aerosols, Atmos. Meas. Tech., 8, 3831–3849, https://doi.org/10.5194/amt838312015, 2015. a
Chance, K. V. and Spurr, R. J. D.: Ring effect studies: Rayleigh scattering, including molecular parameters for rotational Raman scattering, and the Fraunhofer spectrum, Appl. Optics, 36, 5224–5230, 1997. a
Chimot, J., Vlemmix, T., Veefkind, J. P., de Haan, J. F., and Levelt, P. F.: Impact of aerosols on the OMI tropospheric NO_{2} retrievals over industrialized regions: how accurate is the aerosol correction of cloudfree scenes via a simple cloud model?, Atmos. Meas. Tech., 9, 359–382, https://doi.org/10.5194/amt93592016, 2016. a
Chimot, J., Veefkind, J. P., Vlemmix, T., de Haan, J. F., Amiridis, V., Proestakis, E., Marinou, E., and Levelt, P. F.: An exploratory study on the aerosol height retrieval from OMI measurements of the 477 nm O_{2} − O_{2} spectral band using a neural network approach, Atmos. Meas. Tech., 10, 783–809, https://doi.org/10.5194/amt107832017, 2017. a
Cox, A. N. (Ed.): Allen's Astrophysical Quantities, fourth edition, The Cox A. N. book, Springer, New York, NY, https://doi.org/10.1007/9781461211860, 2000. a
De Haan, J. F.: Effect of spectral features in the O_{2} A band on the combined retrieval of cloud properties and NO_{2}, TNTROPOMIKNMI045, 2009. a
De Haan, J. F.: TradeOff between DOAS and Optimal Estimation, RPTROPOMIKNMI051, 2010. a
De Haan, J. F.: DISAMAR: Determining Instrument Specifications and Analyzing Methods for Atmospheric Retrieval – Algorithms and Background, RPTROPOMIKNMI066, 22 January, https://gitlab.com/KNMIOSS/disamar/disamar//blob/main/doc/DISAMAR_V3.1.0_algorithms_20220228.pdf (last access: 5 September 2022), 2012. a, b, c, d, e
De Haan, J. F.: Disamar software v4.1.5 (4.1.5), Zenodo [code], https://doi.org/10.5281/zenodo.6304984, 2022. a
De Haan, J. F., Bosma, P., and Hovenier, J. W.: The adding method for multiple scattering calculations of polarized light, Astron. Astrophys., 181, 371–391, 1987. a, b, c, d, e, f, g
De Rooij, W. A. and van der Stap, C. C. A. H.: Expansion of Mie scattering matrices in generalized spherical functions, Astron. Astrophys., 131, 237–248, 1984. a, b
Frankenberg C., Butz. A., and Toon, G. C.: Disentangling chlorophyll fluorescence from atmospheric scattering effects in O_{2} Aband spectra of reflected sunlight, Geophys. Res. Lett, 38, L03801, https://doi.org/10.1029/2010GL045896, 2011. a, b
Griewank, A. and Walther, A.: Evaluating Derivatives, Principles and Techniques of Algorithmic Differentiation, Society for Industrial and Applied Mathematics, 2nd edition, https://doi.org/10.1137/1.9780898717761, 2008. a
Hansen, J. E. and Travis, L. D.: Light scattering in planetary atmospheres, Space Sci. Rev., 16, 527–610, https://doi.org/10.1007/BF00168069, 1974. a
Hasekamp, O. P. and Landgraf, J.: Linearization of vector radiative transfer with respect to aerosol properties and its use in satellite remote sensing, J. Geophys. Res., 110, D04203, https://doi.org/10.1029/2004JD005260, 2005. a
Hovenier, J. W., van der Mee, C. V. M., and Domke, H.: Transfer of polarized light in planetary atmospheres; basic concepts and practical methods, Kluwer, Dordrecht, 57389048, 2004. a, b, c, d
Joiner, J., Yoshida, Y., Vasilkov, A. P., Yoshida, Y., Corp, L. A., and Middleton, E. M.: First observations of global and seasonal terrestrial chlorophyll fluorescence from space, Biogeosciences, 8, 637–651, https://doi.org/10.5194/bg86372011, 2011. a
Karpowicz, B. M., Stegmann, P. G., Johnson, B. T., Christophersen, H. W., Hyer, E. J., Lambert, A., and Simon, E.: pyCRTM: A python interface for the community radiative transfer model, J. Quant. Spectrosc. Ra., 288, 108263, https://doi.org/10.1016/j.jqsrt.2022.108263., 2022. a
Kroon, M., de Haan, J. F., Veefkind, J. P., Froidevaux, L., Wang, R., Kivi, R., and Hakkarainen, J. J.: Validation of operational ozone profiles from the Ozone Monitoring Instrument, J. Geophys. Res., 116, D18305, https://doi.org/10.1029/2010JD015100, 2011. a, b
Lenoble, J., Herman, M., Deuzé, J. L., Lafrance, B., Santer, R., and Tanré, D.: A successive order of scattering code for solving the vector equation of transfer in the earth’s atmosphere with aerosols, J. Quant. Spectrosc. Ra. 107, 479–507, 2007. a
Lu, C.H., Liu, Q., Wei, S.W., Johnson, B. T., Dang, C., Stegmann, P. G., Grogan, D., Ge, G., Hu, M., and Lueken, M.: The Aerosol Module in the Community Radiative Transfer Model (v2.2 and v2.3): accounting for aerosol transmittance effects on the radiance observation operator, Geosci. Model Dev., 15, 1317–1329, https://doi.org/10.5194/gmd1513172022, 2022. a
Min, Q. and Duan, M.: A successive order of scattering model for solving vector radiative transfer in the atmosphere, J. Quant. Spectrosc. Ra., 87, 243–259, 2004. a
Nanda, S., de Graaf, M., Sneep, M., de Haan, J. F., Stammes, P., Sanders, A. F. J., Tuinder, O., Veefkind, J. P., and Levelt, P. F.: Error sources in the retrieval of aerosol information over bright surfaces from satellite measurements in the oxygen A band, Atmos. Meas. Tech., 11, 161–175, https://doi.org/10.5194/amt111612018, 2018. a, b
Nanda, S., de Graaf, M., Veefkind, J. P., ter Linden, M., Sneep, M., de Haan, J., and Levelt, P. F.: A neural network radiative transfer model approach applied to the Tropospheric Monitoring Instrument aerosol height algorithm, Atmos. Meas. Tech., 12, 6619–6634, https://doi.org/10.5194/amt1266192019, 2019. a, b
Platt, U.: Differential optical absorption spectroscopy (DOAS), in: Air Monitoring by Spectroscopic Techniques...”, Wiley & Sons, New York, USA, 27–84, 1994. a
Platt, U. and Stutz, J.: Differential Optical Absorption Spectroscopy (DOAS) Principles and Applications, Springer, Heidelberg, https://doi.org/10.1007/9783540757764, 2008. a
Rodgers, C. D.: Inverse Methods for Atmospheric Sounding Theory and Practice, World Scientific, Singapore, https://doi.org/10.1142/3171, 2000. a, b, c, d, e
Rothman, L. S., Gordon, I. E., Barbe, A., Benner, D. C., Bernath, P. F., Birk, M., Boudon, V., Brown, L. R., Campargue, A., Champion, J.P., Chance, K., Coudert, L. H., Dana, V., Devi, V. M., Fally, S., Flaud, J.M., Gamache, R. R., Goldman, A., Jacquemart, D., Kleiner, I., Lacome, N., Lafferty, W. J., Mandin, J.Y., Massie, S. T., Mikhailenko, S. N., Miller, C. E., MoazzenAhmadi, N., Naumenko, O. V., Nikitin, A. V., Orphal, J., Perevalov, V. I., Perrin, A., PredoiCross, A., Rinsland, C. P., Rotger, M., Šimec̆ková, M., Smith, M. A. H., Sung, K., Tashkun, S. A., Tennyson, J., Toth, R. A., Vandaele, A. C., and Vander Auwera, J.: The HITRAN 2008 molecular spectroscopic database, J. Quant. Spectrosc. Ra., 110, 533572, https://doi.org/10.1016/j.jqsrt.2009.02.013, 2009. a
Sanders, A. F. J., de Haan, J. F., Sneep, M., Apituley, A., Stammes, P., Vieitez, M. O., Tilstra, L. G., Tuinder, O. N. E., Koning, C. E., and Veefkind, J. P.: Evaluation of the operational Aerosol Layer Height retrieval algorithm for Sentinel5 Precursor: application to O_{2} A band observations from GOME2A, Atmos. Meas. Tech., 8, 4947–4977, https://doi.org/10.5194/amt849472015, 2015. a, b
Saunders, R., Hocking, J., Turner, E., Rayer, P., Rundle, D., Brunel, P., Vidot, J., Roquet, P., Matricardi, M., Geer, A., Bormann, N., and Lupu, C.: An update on the RTTOV fast radiative transfer model (currently at version 12), Geosci. Model Dev., 11, 2717–2737, https://doi.org/10.5194/gmd1127172018, 2018. a
Spurr, R. J. D.: Simultaneous derivation of intensities and weighting functions in a general pseudospherical discrete ordinate radiative transfer treatment, J. Quant. Spectrosc. Ra., 75, 129–175, 2002. a, b
Spurr, R. J. D.: VLIDORT: A linearized pseudospherical vector discrete ordinate radiative transfer code for forward model and retrieval studies in multilayer multiple scattering media, J. Quant. Spectrosc. Ra., 102, 316–342, 2006. a
Stam, D. M., Aben, I., and Helderman, F.: Skylight polarization spectra: Numerical simulation of the Ring effect, J. Geophys. Res., 107, 4419, https://doi.org/10.1029/2001JD000951, 2002. a
Stammes, P., de Haan, J. F., and Hovenier, J. W.: The polarized internal radiation field of a planetary atmosphere, Astron. Astrophys., 225, 239–259, 1989. a, b
Stegmann, P. G., Johnson, B. T., Moradi, I., Karpowicz, B., and McCarty, W.: A deep learning approach to fast radiative transfer, J. Quant. Spectrosc. Ra., 280, https://doi.org/10.1016/j.jqsrt.2022.108088, 2022. a
Tran, H. and Hartmann, J.M.: An improved O_{2} A band absorption model and its consequences for retrievals of photon paths and surface pressures, J. Geophys. Res., 113, D18104, https://doi.org/10.1029/2008JD010011, 2008. a
Van de Hulst, H. C.: Light scattering by small particles, Wiley (also Dover), https://doi.org/10.1002/actp.1984.010350426, 1981. a
Vasilkov, A. P., Joiner, J., Gleason, J., and Bhartia, P. K.: Ocean Raman scattering in satellite backscatter UV measurements, Geophys. Res. Lett., 29, 1837, https://doi.org/10.1029/2002GL014955, 2002. a
Veefkind, J. P., Aben, I., McMullan, K., Förster, H., de Vries, J., Otter, G., Claas, J., Eskes, H. J., de Haan, J. F., Kleipool, Q., van Weele, M., Hasekamp, O., Hoogeveen, R., Landgraf, J., Snel, R., Tol, P., Ingmann, P., Voors, R., Kruizinga, B., Vink, R., Visser, H., and Levelt, P. F.: TROPOMI on the ESA Sentinel5 Precursor: A GMES mission for global observations of the atmospheric composition for climate, air quality and ozone layer applications, Remote Sens. Environ., 120, 70–83, https://doi.org/10.1016/j.rse.2011.09.027, 2012. a
Veefkind, J. P., Keppens, A., and de Haan, J. F.: TROPOMI Ozone Profile ATBD, Report S5PKNMIL20004RP, KNMI, 20211022, https://sentinel.esa.int/documents/247904/2476257/Sentinel5PTROPOMIATBDOzoneProfile.pdf (last access: 2 September 2022), 2021. a
 Abstract
 Introduction
 Novel features of DISAMAR
 DISAMAR forward model
 DISAMAR retrieval methods
 Examples of applications
 Summary and concluding remarks
 Appendix A: Derivation of the derivative formulae by means of linearization
 Appendix B: Pseudospherical correction
 Appendix C: Conversion of pressure grid to altitude grid
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Novel features of DISAMAR
 DISAMAR forward model
 DISAMAR retrieval methods
 Examples of applications
 Summary and concluding remarks
 Appendix A: Derivation of the derivative formulae by means of linearization
 Appendix B: Pseudospherical correction
 Appendix C: Conversion of pressure grid to altitude grid
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References