HAL Id: insu-02183026
https://hal-insu.archives-ouvertes.fr/insu-02183026
Submitted on 8 Dec 2020
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Improved ozone DIAL retrievals in the upper
troposphere and lower stratosphere using an optimal
estimation method
Ghazal Farhani, Robert J. Sica, Sophie Godin-Beekmann, Gérard Ancellet,
Alexander Haefele
To cite this version:
Ghazal Farhani, Robert J. Sica, Sophie Godin-Beekmann, Gérard Ancellet, Alexander Haefele.
Improved ozone DIAL retrievals in the upper troposphere and lower stratosphere using an opti-
mal estimation method. Applied optics, Optical Society of America, 2019, 58 (6), pp.1374-1385.
�10.1364/AO.58.001374�. �insu-02183026�
Research Article Applied Optics 1
Improved ozone DIAL retrievals in the upper
troposphere and lower stratosphere using an optimal
estimation method
GHAZAL FARHANI
1,*
, ROBERT J. SICA
1
, SOPHIE GODIN-BEEKMANN
2
, GÈRARD ANCELLET
2
, AND
ALEXANDER HAEFELE
3
1
Department of Physics and Astronomy, The University of Western Ontario, 1151 Richmond St., London, ON, N6A 3K7
2
Observatoire de Versailles Saint-Quentin-en-Yvelines, Guyancourt, France 78280
3
Federal Office of Meteorology and Climatology MeteoSwiss, Payerne, Switzerland 1530
*
Corresponding author: gfarhani@uwo.ca
Compiled October 31, 2018
We have implemented a first-principle Optimal Estimation Method to retrieve ozone density profiles us-
ing simultaneously tropospheric and stratospheric Differential Absorption Lidar (DIAL) measurements.
Our retrieval extends from 2.5 km to about 42 km in altitude, and in the upper troposphere and the lower
stratosphere (UTLS) it shows a significant improvement in the overlapping region, where the OEM can
retrieve a single ozone profile consistent with the measurements from both lidars. Here stratospheric
and the tropospheric measurements from the Observatoire de Haute Provence are used, and the OEM
retrievals in the UTLS region compared with coincident ozonesonde measurements. The retrieved ozone
profile have a small statistical uncertainty in the UTLS region relative to individual determinations of
ozone from each lidar, and the maximum statistical uncertainty does not exceed a maximum of 7%. © 2018
Optical Society of America
http://dx.doi.org/10.1364/ao.XX.XXXXXX
1. INTRODUCTION
The upper troposphere and lower stratosphere (UTLS) extends from about 6 km to 25 km in height and plays
a significant role in the atmospheric climate system. In this region of the atmosphere, even small changes
in temperature and in the distribution and concentration of greenhouse gases can result in large changes in
atmospheric radiative forcing, which can trigger climate change [1, 2].
Ozone in the upper troposphere acts as the third largest greenhouse gas contributing to the radiative
forcing of climate change [
1
,
3
]. The ozone distribution in the UTLS is the result of transport mechanisms and
photochemical reactions. Because of stratospheric tropospheric exchange, large spatial and temporal variability
can be observed in the UTLS [4].
In many studies on the UTLS ozone, satellite-borne instruments are used. In limb-viewing instruments,
the elevation angle of the line-of-sight varies during the measurements. As a result, limb sounders can
provide good vertical resolution (about 2 km to 4 km). However, at lower altitudes (lower troposphere), the
atmosphere becomes nearly opaque, and the limb-viewing instruments have difficulties measuring trace
gases. Nadir-viewing instruments can provide measurements in the lower troposphere, but their vertical
resolution is limited (about 6 km to 7 km). Occultation instruments use the Sun or other stars as the source
of radiation, and they can obtain measurements with higher vertical resolution (about 1 km to 2 km). Solar
occultation instruments are restricted by the number of sunsets and sunrises they encounter in one orbit, while
Research Article Applied Optics 2
stellar occultation instruments are limited by the weakness of the stellar source compared to the Sun. The
combination of measurements from different geometrical-based satellite instruments has been used to measure
ozone density.
In the UTLS, large biases (differences) in ozone measurements are reported between instruments. Although
the difference between the data sets is more significant in the tropics and high latitudes (about
±
30%),
significant bias exists at mid-latitudes (about
±
10%) [
5
]. Therefore, a continued detailed intercomparison
between satellite instruments, as well as between satellites and other instruments, is needed, including both
airborne measurements and ground-based measurements.
Differential Absorption Lidar (DIAL) systems provide ozone measurements with high vertical and temporal
resolutions. For example, observatories such as the Canadian Network for the Detection of Atmospheric
Change (CANDAC) Polar Environment Atmospheric Research Laboratory (PEARL) in Eureka, Maïdo ob-
servatory in Reunion Island, the Observatoire de Haute Provence (OHP) in France, and the NASA Table
Mountain Observatory (TMO) in the United States are equipped with both tropospheric and stratospheric
lidars. At the Eureka observatory, the tropospheric lidar system makes measurements from 0.5 km to about
8 km in altitude and the stratospheric lidar system operates from about 4 km to 35 km [
6
,
7
]. At the Maïdo
observatory, the tropopspheric DIAL makes measurements from 6 km to 16 km, and the stratospheric DIAL
operates in the 13 km to 38 km region [
8
]. At the OHP observatory, the tropospheric DIAL system operates
from 2.5 km to about 14.5 km, and the stratospheric DIAL operates from about 10 km to 45 km [
9
,
10
]. At the
TMO, the tropospheric DIAL system obtains measurements from 3 km to 18 km, and the stratospheric DIAL
system from 10 km to 40 km [
11
,
12
]. Although these systems can produce satisfactory ozone profiles in their
overlapping region (from tropospheric lidar to stratospheric lidar), the uncertainty of merging is not well
defined. Providing a single ozone profile with a full uncertainty budget using both sets of measurements can
significantly improve our measurements of ozone in the UTLS [13–15].
Here we apply the Optimal Estimation Method (OEM) to tropospheric and stratospheric DIAL measure-
ments. Measurements from these two systems are simultaneously used by the retrieval to obtain a single
ozone profile. Using the OEM there is no need to “merge” or “glue” level 0 profiles. Moreover, the input
measurements can be in different units with different measurement grids (for example a mix of analog and
digital measurements). Additionally, a full uncertainty budget, including both the systematic and statistical
uncertainties, is calculated for each individual profile. The OEM also provides averaging kernels of the re-
trievals, which allows comparison of the profiles with other measurements which can account for differences in
vertical resolution, such as when compared to space-based measurements. Other atmospheric and systematic
parameters such as air density, the dead time of the system, and the background counts can be retrieved along
with ozone profiles. The application of OEMs to aerosol lidar measurements, Rayleigh scatter temperatures,
and Raman scatter water vapour retrievals has been studied and discussed in detail [
16
–
18
]. In addition, we
have recently demonstrated an OEM for DIAL stratospheric ozone retrievals [
19
], which we will now expand
to include measurements from tropospheric ozone DIAL systems.
In this paper, focusing on the UTLS region, we show a first principle OEM to retrieve a single ozone
profile by using both tropospheric and stratospheric DIAL measurements directly from the raw (level 0)
measurements using the lidar equation as the forward model. In Section 2, pre-processing steps prior to
applying the traditional DIAL algorithm, as well as the OEM, are discussed. Moreover, the state vectors and
the
b
parameter quantities are defined and a brief overview of the lidar’s specifications is given. In Section 3
results of the OEM retrieval, using both tropospheric and stratospheric lidar measurements, are discussed
in detail. In this Section, we also show our results of comparison between the ozonesonde profiles and our
retrievals. Section 4 is the summary of the paper, and in Section 5 we discuss our future plans. Details of how
to apply our method to a standalone tropospheric DIAL are given in the Appendix.
2. METHODOLOGY
In the DIAL system, two wavelengths are simultaneously transmitted to the atmosphere. One of the emitted
wavelengths is strongly absorbed by the constituent of interest (called the “on-line” wavelength) and the other
is weakly absorbed (called the “off-line” wavelength). For ozone measurements, selecting a wavelength pair
Research Article Applied Optics 3
depends on the altitude range of the measurements. For most studies, the ultraviolet (UV) spectrum is the
most efficient spectral region. A pair of wavelengths with strong UV absorption is needed to detect the small
amount of ozone which resides in the troposphere. However, for stratospheric ozone measurements, choosing
a laser that can reach higher altitudes in the stratosphere is the main concern [11, 20, 21].
The traditional analysis method for ozone number density uses the derivative of the ratio between the
“on-line” and “off-line” channels to calculate the ozone number density
n
o
3
(z)
. A detailed discussion on the
tropospheric and stratospheric ozone retrievals can be found in [
10
,
22
–
25
]. In the traditional analysis, some
corrections are applied to the raw lidar measurements, for example background counts should be removed.
In many systems this requires including the effects of signal-induced-noise (SIN). Any corrections due to
nonlinearity of the counting system (because of saturation) should also be applied to the raw counts. Finally,
the signals from different channels need to be merged to form a single measurement profile. This “corrected”
count profile is then used to calculate the ozone number density profiles. With the OEM, a forward model
encapsulates the geophysical properties and instrumental characteristics of the system, and our OEM retrieval
uses the raw (level 0) measurements from all available channels. A comprehensive explanation of the OEM
can be found in [26]; a brief description of the OEM follows below.
In the OEM a forward model is defined as the relation between the measurements vector
y
= (
y
1
, y
2
, ..., y
n
),
and the state vector x = (x
1
, x
2
, ..., x
n
). The forward model is:
y = F(x, b) + e (1)
where
b
are the forward model parameters, which are assumed to be known, and
e
is the measurement noise.
We use the lidar equation as the forward model, where the raw counts are the measurements. The lidar
equation for unsaturated counts, N
true
, is:
N
obs
(z, λ
i
) =
C(λ
i
)O(z)
z
2
β(λ
i
, z) exp[−2
Z
∞
0
[σ
O
3
(λ, T(z))n
O
3
(z) + α(λ, z) +
∑
e
σ
e
(λ)n
e
(z)]dz] + N
b
(z, λ
i
) (2)
where
N
obs
(
z, λ
i
) is the number of backscattered photons.
C
(
λ
i
) is the lidar constant, which contains the area
of the receiving telescope, the total efficiency of the lidar system, and energy of the scattered photon. The
geometrical overlap is
O
(
z
), and
β
(
λ
i
, z
) are the atmospheric backscattering coefficients which includes both
molecular and aerosol terms. The first term inside the integral corresponds to ozone absorption in which
σ
O
3
(
T
(
z
)
, λ
i
) is the ozone absorption cross section, which depends on atmospheric temperature, and
n
O
3
(
z
) is
the ozone number density. The second term of the integral,
α
(
λ, z
) contains the extinction coefficient which
is the sum of the extinction due to molecules and particles, and the last term
∑
e
σ
e
(
λ
)
n
e
(
z
) is the extinction
by other absorbers. For ozone studies, the most common interfering gases are
SO
2
,
NO
2
and O
2
. The effect
of O
2
is only considered when the selected “on-line” laser wavelength is shorter than 294 nm [
12
,
27
]. In the
case of heavy volcanic eruption,
SO
2
and
NO
2
can significantly affect ozone retrievals [
28
]. However, in most
cases, for both stratospheric and tropospheric ozone studies the effect of these gases in final ozone retrievals is
negligible. Thus, the last term of integration is typically neglected [10].
The background counts are written as
N
b
(
z
). In the presence of SIN, the background is fitted to an exponential
function of the form:
N
b
(z) = a exp(−bz) + c (3)
where
a, b,
and
c
are coefficients of the fit, which in the traditional method are determined analytically, but are
retrieved in our OEM retrieval using the analytic values as a priori coefficients [29].
When the intensity of the backscattered signal is high, the counting system can be affected by saturation.
This saturation can result in an observed count rate which is less than the true count rate. For a paralyzable
detector, true counts are related to the observed counts N
obs
as follows:
N
obs
= N
true
exp(−κN
true
) (4)
and, for non-paralyzable detectors, the following equation can be used:
N
obs
=
N
true
1 + κN
true
(5)
Research Article Applied Optics 4
where
κ
is the dead time of the detecting system. For the OEM retrieval the value of the dead time for each
channel is retrieved.
A. Implementing the OEM for the OHP lidars
Knowing the measurements vector and its covariance matrix
S
e
, and using an a priori profile and its associated
covariance matrix S
a
, the OEM calculates an optimal a posteriori state by minimizing a cost function:
Cost = (y − Kx)
T
S
−1
y
(y − Kx) + (x −
b
x)
T
S
−1
a
(x −
b
x) (6)
As our forward model is nonlinear, an iterative numerical method is used. For our problem the Levenberg-
Marquardt iteration is a suitable numerical method. Then, the optimized state vector x is given as:
x
i+1
= x
i
+ [(1 + γ
i
)S
−1
a
+ K
T
i
S
y
K
T
i
]
−1
(
[K
T
i
S
−1
y
(y − F(x
i
, b)] − S
−1
a
(x
i
− x
a
)
)
(7)
where
i
is the iteration term,
x
a
is the a priori profile, and
K
=
dF
dx
is the linearisation term for our nonlinear
forward model, called the Jacobian matrix. Finally,
γ
i
is a damping factor for the iteration, which is chosen at
each step to minimize the cost function. As suggested by [
30
] if the value of the cost function increases in a
step,
γ
i
will increase by a factor of 10, and if the value of the cost function decreases in a step,
γ
i
will decrease
by a factor of 2. The iteration stops when the cost function decreases to a value much smaller than the number
of measurements. There are other criteria which result in ceasing the iteration. Further details can be found in
[26].
To understand how measurements and a priori profiles contribute in the final retrievals, an averaging kernel
can be used. The relation between the retrieved state and the true state is described by the averaging kernel of
the retrieval. The averaging kernel is calculated as:
A =
d
b
x
dx
= [K
T
S
−1
y
K + S
−1
a
]
−1
K
T
S
−1
y
K (8)
The retrieved quantity (
b
x) can be written as follows:
b
x = (I − A)x
a
+ Ax + e
r
(9)
where
e
r
is the retrieval uncertainty and
I
is a unity matrix. A perfect retrieval, in the sense all the information
comes from the measurement with no effect from the a priori state, has averaging kernels equal to one, where
the first term of the above equation becomes zero. The width of the averaging kernel gives the resolution of
the retrieval at each height, here defined as the Full Width Half Maximum (FWHM) of each averaging kernel.
In order to find the state vector (from Eq. 7) the following quantities should be known: the measurements
and their covariances, the a priori profiles, the a priori profile’s covariance, and the model (
b
) parameters. The
b
parameters are quantities in the forward model that are not being retrieved, because they are either well-known
or retrieving them is not possible. The uncertainty associated with the retrieval due to the
b
parameters is
calculated after the last iteration of the solution. The forward model and the Jacobians (
K
) for each of the state
vectors are calculated, and the Qpack package is used to perform the retrieval. Details of the Qpack software
are given in [31].
Here we retrieve the ozone density profile, relative air density, dead time values, and background counts.
Overlap functions, ozone cross sections, and Rayleigh scattering cross sections are considered as
b
parameters
in the forward model. Below, we discuss our choices of a priori profiles and
b
parameter values. The covariance
matrices associated with the measurements and a priori profiles are discussed as well, and these values are
summarized in Table 1.
In photon counting mode, when the signal is linear, the measurements statistical uncertainty follows a
Poisson distribution, and the number of counts at each altitude represents the measurement’s variance at
that height. There is no correlation between the digital counts in different layers of the atmosphere, so
the off-diagonal elements of the measurement’s covariance matrix are zero. However, for the OHP lidars,
analog measurements do not follow Poisson distributions. Calculating the measurement variance for each
