# Verification for Different Contrail Parameterizations Based on Integrated Satellite Observation and ECMWF Reanalysis Data

Authors: Jinglin Zhang; Jian Shang; Guoyu Zhang
Journal: Advances in Meteorology
 (2017)

---

## Abstract

Aviation induced cloud termed contrail plays a more and more important role in the climate change, which makes a significant contribution to anthropogenic climate forcing through impacting the coverage of cirrus in the intersection of troposphere and stratosphere. In this paper, we propose one novel automatic contrail detecting method based on Himawari-8 stationary satellite imagery and two kinds of potential contrail coverage (PCC1 and PCC2) from contrail parameterization in ECHAM4 and HadGEM2. In addition, we propose one new climatological index called contrail occurrence and persistence (COP). According to the algorithm identification (AI) and artificial visual inspection (AVI), COP measured from Himawari-8 stationary satellite imagery is related to upper tropospheric relative humidity over ice (RHI) computed with the ECMWF reanalysis data by simple linear regression. Similarly, we compared the linear correlation between COP and PCCs fractions and found that PCC1 has better correspondence with COP than PCC2.


---

## Body

## 1. Introduction

The cloud is an important regulator of the energy budget of the earth system and the most significant uncertainty factor in climate change [1–4]. Contrail is a special kind of cloud, caused by the aircraft engine emissions from the exhaust gas mixed with the surrounding air and the condensation of water vapor. As one kind of condensation trails, contrail is the direct evidence for the impact of air traffic on cloud cover. The ambient relative humidity is close to saturation, and the atmospheric state is relatively stable to produce contrail. In general, the formation of cirrus requires temperature below −40°C and relative humidity of 145%–165% or more in the stable atmospheric state [5], while the contrail can persist in weakly saturated condition (relative humidity 100%–110%). Compared to the natural cirrus, the contrails are more likely to occur and persist at the upper troposphere and lower stratosphere, which affect the coverage of the high-level cirrus. Similar to natural cirrus clouds, it consists of tiny ice particles and interacts with the radiation field, which has a net warming effect on the earth-atmosphere system [6–8]. According to the recent study [9], aircraft induced contrails contribute more than the CO2 caused by the airplane flight to global warming. Thus, one considerate detection framework will help researchers to understand the influence of contrails in climate change better. Generally, contrails are divided into three types: young, mature, and old in terms of contrail width and shape. A young contrail was defined as a linear shape with approximately a width less than 3 km and a maximum distance less than 200 km. The mature contrail maintains a linear shape with an approximate width less than 7 km. The old contrail gradually loses the linear shape and develops into the contrail cirrus. It can be observed that a mature contrail has larger brightness temperature difference (BTD) values than old contrails.For the human visual verification, it is easy to identify linear structures from both ground-based and satellite images as shown in Figure1. Thus, early explorations about contrail detection mainly exploit the visual examination to detect contrails in satellite images. For example, Degrand et al. (1990) applied hardcopy displays of high-resolution Defense Meteorological Satellite Program data to identify contrails. Degrand et al. (1991) manually interpreted contrails on satellite images over North America and in particular over the United States. Travis (1996) manually counted pixels to determine the width and length of the contrails. All of these manual inspection methods are highly subjective and very time-consuming.Figure 1
Contrails ground-based and satellite observation.
(a)
Ground-based observation of contrails
(b)
Himawari-8 satellite BTD image (band(13)–band (15)) with 2 km spatial resolution. 2017/4/5 UTC-03:50, 115°E–144°E, 30°N–44°NIn order to overcome the limitations of manual detection methods for contrails, researchers attempted to develop automatic algorithms to detect contrails depending on satellite imagery. Mannstein et al. (1999) introduced a contrail detection algorithm (CDA) to recognize linear-shaped contrails on satellite images, which is operated on the brightness temperature difference between the 10.8μm and 12 μm channels. Recently, other studies incorporated CDA to detect contrails. Palikonda et al. (2001) identified the contrails over the United States in AVHRR and MODIS images with CDA method. More recently, Hetzheim [10] proposed a very complicated approach using mathematical methods to detect contrails based on texture or stochastic behaviors of contrails. Although these mathematical methods may distinguish contrails well from background of the satellite images, they are very time-consuming.To the best of our knowledge, there are still some open questions concerning contrail occurrence, persistence, and spreading, which cause uncertainties in the contrail detection. New generation meteorological satellite like Himawari-8 with 16 multispectral channels captures high quality satellite imagery. Meanwhile stationary satellite imagery with both high spatial and temporal resolution are necessary to locate the linear contrails and track the process from linear contrail (occurrence) to contrail cirrus (persistence), in order to resolve the problem of contrail cirrus overlap and limitation of polar satellite overpass time relative to traffic diurnal cycle. With automatic contrail detection methods, previous research works [11–13] mainly focused on heavily flown regions like the North Atlantic flight corridor and Europe regions to compute the contrail coverage and related radiative forcing impact. Little study about contrail was carried out in the Asian region, especially the east coast of China and Japan.Besides the satellite observation, there also have been many studies about the simulation of contrail formation and persistence through the parameterization in a general circulation model (GCM). Early, Schmidt [14] and Appleman [15] provided the thermodynamic explanation of the contrail formation process. According to the Schmidt-Appleman criterion (SAC), further studies [16, 17] pointed out that contrails are usually short-lived once they form in dry ambient air, and if the ambient relative humidity with respect to ice exceeds ice saturation, they can persist at few hours and spread with the ambient wind field and may develop into cirrus clouds. Afterward, Sausen et al. [18] diagnosed the global seasonal contrails coverage from assimilated meteorological data according to a method in a GCM. This method can calculate high cloud coverage from ambient temperature and relative humidity and also can apply SAC to judge the contrail existence when there are aviation activities. A. Schmitt and B. Brunner had provided the data for this method; they use the European Centre for Medium-Range Weather Forecasts (ECMWF) as the meteorological data and the DLR-2 (Deutschen Zentrums für Luft- und Raumfahrt, DLR) data set [19] as the aviation data. The following study Ponater et al. [20] presented the simulated contrail properties and the radiative forcing distribution resulting from contrails. Marquart et al. [21] continually investigated the future development of contrails within a consistent framework, including contrail coverage, contrail properties, and radiative forcing but still limited to line-shaped contrails. Nowadays, the simulation of contrails in a GCM has developed into containing line-shaped contrails and contrail cirrus, while the simulation of contrail coverage or contrail cirrus coverage depends on the calculation of potential contrail coverage (PCC, [22]). The PCC is a method applied to determine the ability of the ambient air to allow contrail formation and consists of SAC and the parameterization of natural cirrus coverage from temperature and humidity in a GCM. Burkhardt and Kärcher [9, 23] and Burkhardt et al., 2015, developed a process-based contrail cirrus module in ECHAM4 to determine contrail formation, persistence, and translation and has used their module to simulate contrail cirrus coverage, the radiative forcing, and related changes in the natural cirrus clouds. Similarly, Rap et al. [24] developed a contrail parameterization by adapting the contrail parameterization of Ponater et al. [20] to the second version of the UK Hadley Centre Global Environmental Model (HadGEM2).Therefore, we propose one novel object-based contrail detecting method. Based on high temporal resolution imagery of Himawari-8 stationary satellite, it not only detects the contrail occurrence but also reflects the process of contrails persistence and spreading. Different from previous works, we define one new climatological index termed contrail occurrence and persistence (COP). In our statistics experiments, we calculate the fraction of COP with algorithm identification (AI) and artificial visual inspection (AVI), respectively, while analyzing the linear correlation between COP fraction, relative humidity over ice (RHI), and PCC. This paper is organized as follows: Section2 presents the Himawari-8 satellite and reanalysis data we used. Section 3 illustrates the proposed contrail detection method, the way to calculate the COP index, and PCC parameterization. Section 4 introduces our experimental result and discussion. A brief conclusion is given in Section 5.


## 2. Data

### 2.1. Himawari-8 Satellite

Himawari-8 stationary satellite was launched in Oct. 2014. The Advanced Himawari Imager (AHI) on board Himawari-8 is a 16-channel multispectral imagery to capture visible light and infrared images of the Asia-Pacific region. The instrument was designed and built by Exelis Geospatial Systems and has similar spectral and spatial characteristics to the Advanced Baseline Imager (ABI) planned for use in the American GOES-R satellites. The AHI produces images with a resolution down to 500 m and provides full disk observations every 10 mins and the whole Japan territory every 2.5 minutes. Table1 introduces the characteristics of the Himawari-8 AHI instrument with different channels.Table 1
Characteristics of the Himawari-8 AHI instrument.
Band
Center wavelength
Bandwidth
Resolution at SSP(subsatellite point)
Prime measurement objectives and use of sample data
(
1
)
470 nm
50 nm
1.0 km
Daytime aerosol over land, coastal water mapping
(
2
)
510 nm
20 nm
1.0 km
Green band to produce color composite imagery
(
3
)
644 nm
30 nm
0.5 km
Daytime vegetation/burn scar and aerosols over water, winds
(
4
)
860 nm
20 nm
1.0 km
Daytime cirrus cloud
(
5
)
1610 nm
20 nm
2.0 km
Daytime cloud-top phase and particle size, snow
(
6
)
2260 nm
20 nm
2.0 km
Daytime cloud properties, particle size, vegetation, snow
(
7
)
3.90µm
0.22µm
2.0 km
Surface and cloud, fog at night, fire, winds
(
8
)
6.185µm
0.37µm
2.0 km
High-level atmospheric water vapor, winds, rainfall
(
9
)
6.95µm
0.12µm
2.0 km
Mid-level atmospheric water vapor, winds, rainfall
(
10
)
7.40µm
0.17µm
2.0 km
Lower-level water vapor, winds and SO2
(
11
)
8.50µm
0.32µm
2.0 km
Total water for stability, cloud phase, dust, SO2, rainfall
(
12
)
9.61µm
0.18µm
2.0 km
Total ozone, turbulence, winds
(
13
)
10.35µm
0.30µm
2.0 km
Surface and cloud
(
14
)
11.20µm
0.20µm
2.0 km
Imagery, SST, clouds, rainfall
(
15
)
12.30µm
0.30µm
2.0 km
Total water, ash, SST
(
16
)
13.30µm
0.20µm
2.0 km
Air temperature, cloud heights and amounts


### 2.2. Reanalysis Data

In this paper, we make use of the ECMWF reanalysis data (ERA-Interim, spring in 2016, D. P. Dee et al., 2011). Monthly values of temperature, specific humidity, and ice water content are employed. Horizontally, we have chosen the resolution of variables to a0.25°×0.25° latitude-longitude grid, and vertically we choose 200, 225, 250, and 300 hPa levels, those included the most altitudes of aircraft levels.


## 2.1. Himawari-8 Satellite

Himawari-8 stationary satellite was launched in Oct. 2014. The Advanced Himawari Imager (AHI) on board Himawari-8 is a 16-channel multispectral imagery to capture visible light and infrared images of the Asia-Pacific region. The instrument was designed and built by Exelis Geospatial Systems and has similar spectral and spatial characteristics to the Advanced Baseline Imager (ABI) planned for use in the American GOES-R satellites. The AHI produces images with a resolution down to 500 m and provides full disk observations every 10 mins and the whole Japan territory every 2.5 minutes. Table1 introduces the characteristics of the Himawari-8 AHI instrument with different channels.Table 1
Characteristics of the Himawari-8 AHI instrument.
Band
Center wavelength
Bandwidth
Resolution at SSP(subsatellite point)
Prime measurement objectives and use of sample data
(
1
)
470 nm
50 nm
1.0 km
Daytime aerosol over land, coastal water mapping
(
2
)
510 nm
20 nm
1.0 km
Green band to produce color composite imagery
(
3
)
644 nm
30 nm
0.5 km
Daytime vegetation/burn scar and aerosols over water, winds
(
4
)
860 nm
20 nm
1.0 km
Daytime cirrus cloud
(
5
)
1610 nm
20 nm
2.0 km
Daytime cloud-top phase and particle size, snow
(
6
)
2260 nm
20 nm
2.0 km
Daytime cloud properties, particle size, vegetation, snow
(
7
)
3.90µm
0.22µm
2.0 km
Surface and cloud, fog at night, fire, winds
(
8
)
6.185µm
0.37µm
2.0 km
High-level atmospheric water vapor, winds, rainfall
(
9
)
6.95µm
0.12µm
2.0 km
Mid-level atmospheric water vapor, winds, rainfall
(
10
)
7.40µm
0.17µm
2.0 km
Lower-level water vapor, winds and SO2
(
11
)
8.50µm
0.32µm
2.0 km
Total water for stability, cloud phase, dust, SO2, rainfall
(
12
)
9.61µm
0.18µm
2.0 km
Total ozone, turbulence, winds
(
13
)
10.35µm
0.30µm
2.0 km
Surface and cloud
(
14
)
11.20µm
0.20µm
2.0 km
Imagery, SST, clouds, rainfall
(
15
)
12.30µm
0.30µm
2.0 km
Total water, ash, SST
(
16
)
13.30µm
0.20µm
2.0 km
Air temperature, cloud heights and amounts


## 2.2. Reanalysis Data

In this paper, we make use of the ECMWF reanalysis data (ERA-Interim, spring in 2016, D. P. Dee et al., 2011). Monthly values of temperature, specific humidity, and ice water content are employed. Horizontally, we have chosen the resolution of variables to a0.25°×0.25° latitude-longitude grid, and vertically we choose 200, 225, 250, and 300 hPa levels, those included the most altitudes of aircraft levels.


## 3. Method

As shown in Figure2, our research region is sorted into 1°×1° degree grid boxes extending from 115°E to 144°E and 30°N to 44°N. Three scenes of the Himawari-8 images with different time scales are illustrated in Figure 3, which represent different scenarios to evaluate the proposed contrail detection method. Table 2 presents the information of these Himawari-8 scenes: date and time, spatial extent, and a brief description of contrails.Table 2
Different scenarios of contrails distribution.
ID
Date and time
Spatial extent
Scenario description
(a)
2017/4/9, 03:20/03:30/03:40
118°E–120°E, 34°N–36°N
Contrails cross with each other
(b)
2017/4/4, 03:30/03:40/03:50
128°E–130°E, 42°N–44°N
Contrails lie closely with each other
(c)
2017/3/31, 00:10/00:40/01:00
126°E–132°E, 32°N–38°N
Contrails are densely distributed in the imageFigure 2
BTD image (band(13)–band (15)) with 1°×1° degree grid boxes (115°E–144°E, 30°N–44°N).Figure 3
Different scenarios (a), (b), and (c) of contrails distribution according to Table2.
(a)
(b)
(c)


### 3.1. Contrail Detection Method

Different from traditional CDA method and its variants only considering the spectral information of the image, our work focuses on the object-based method for the detecting contrails in each grid box of Himawari-8 images as shown in Figure2. The main process consists of the following two steps.


#### 3.1.1. Multiresolution Segmentation

In our method, objects in Himawari-8 images are divided into different groups based on contrast, similarity, and shape of object, which are termed heterogeneity. There are several basic steps in the multiresolution segmentation procedure.(1)
The first step of segmentation procedure is to define a single pixel as the starting point.(2)
The single pixel fuses with other pixels as one group, if they are within the threshold of heterogeneity criteria.(3)
Pixel groups fuse with those neighbor groups.(4)
All of the pixels have been assigned to different objects until no further fusing is found.


#### 3.1.2. Object-Based Classification

A nearest neighbor classification method is adopted to identify the Himawari-8 image into two classes: contrail and noncontrail object after multiresolution segmentation:(1)
It assigns several samples with contrail and noncontrail label, respectively. Objects should be verified with shape, size, and spectral characteristics.(2)
Contrail and noncontrail samples are trained from the segmented objects. Contrail and noncontrail classes are linked with the segmented image objects using the k-nearest neighbor classification method.(3)
Classification returns a fuzzy membership value for all the segmented image objects. A value of 0 means the object does not belong to contrail class, while a value of 1 means the object most likely belongs to contrail class.


### 3.2. Artificial Visual Analysis

With the high temporal resolution stationary satellite images, we not only explore the occurrence of linear contrail but also study the development from linear contrail to contrail cirrus. In our evaluation, we investigate the satellite images of spring 2016 (March, May, and April) in our statistics analysis. After a large amount of observation analysis, we choose 30 min as the time scale of statistic, so as to avoid the overestimation of contrail, as well as fully considering the spreading of contrails. Thus, 4416 satellite images (24 per day and 92 days) are analyzed with both the proposed contrail detection method and artificial visual inspection. In the previous studies [1, 25–27], the reconstruction of contrail occurrence fraction adopts commercial flight track data to mask out regions of no air traffic. Their work has obvious limitations as follows:(1)
The polar satellite overpasses the observed region twice per day, which cause the severe insufficiency in observation sample.(2)
This kind of reconstruction did not consider the contrail spreading and persistence, with limited observation sample.Based on a very large amount of observation, some contrails could persist more than 5 h and spread 400–500 km at most, according to the annual mean wind speed at 200–300 hPa. Thus we define one new climatological index, termed contrail occurrence and persistence (COP). First, each of the grid boxes in satellite images is divided into 8 × 8 subgrid boxes; the COP is the ratio between the number of subgrids that contain detected contrails and the total number of subgrids, as illustrated in Figure4. Different from the traditional climatological index of contrail based on limited satellite overpass images, the COP index will count and aggregate detected contrails with 30 min time interval, which will reflect the occurrence and persistence of contrail.Figure 4
Illustration of8×8 subgrid box with detected linear contrails. The contrails occurrence fraction defined in this example is (8 + 8 + 4 + 3)/64 = 0.36.


### 3.3. Relative Humidity over Ice (RHI)

The contrails occurrence fraction presents a significant linear correlation versus the RHI on North America [25]. Therefore, for the comparison between the contrails occurrence fraction versus RHI on East Asia, we use the reanalysis data of temperature to calculate RHI by a formula as follows:(1)RHI=eesatiT≈qqsati,where q is the specific humidity, qsat(i) is the saturation specific humidity over ice, and esatiT is the saturation pressure of water vapor over ice and calculated by an empirical formula as follows:(2)esatiT=6.108exp⁡21.875T-276.16T-7.66,


### 3.4. Potential Contrail Coverage (PCC) Parameterization

Contrail formation depends on the surrounding temperature, moisture, and pressure fields that create an enough cold and moist ambient air for the mixing between the aircraft exhaust plume and its environment. Besides, we use SAC to diagnose the condition of ambient air with persistent contrails [16, 17]. The PCC is a method applied to determine the ability of the ambient air to allow contrail formation. An introduction necessary for PCC to estimate contrail coverage requires the fractional area of the atmosphere to support a persistent contrail formation [20]. Hence, in this paper, we adopt a similar method of PCC parameterization developed by Ponater et al. [20] for the ECHAM4 climate model. Relying on the thermodynamic theory, an empirical formula of the threshold values of temperature (Tcontr) and relative humidity (rcontr(T)) supporting contrail formation and persistence has been given according to Schumann [16]. In addition, contrail formation and persistence depend on some specific parameters and variables including the ambient temperature (T), pressure (p), specific humidity (q), the emission index of water vapor (EIH2o), the specific heat of fuel combustion (Q), and the propulsion efficiency of the aircraft (η). The thermodynamic theory predicts threshold values for temperature (Tcontr) and humidity (rcontrT). No contrails are possible if T>Tcontr and r<rcontrT. The formula of Tcontr (in °C) is given by the following:(3)Tcontr=-46.46+9.43ln⁡G-0.053+0.72ln⁡G-0.0532,where G is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute temperature versus water vapor partial pressure diagram and has the unit of Pa K−1, defined as follows:(4)G=EIH2o·cpmpεQ1-η,where EIH2o=1.25, ε=0.622 is the ratio of molecular masses of water and dry air, Q=4.3×107 J kg−1 K−1, η=0.3, and cpm is the isobaric heat capacity of moist air, which has the unit of J kg−1 K−1 and is given by the following:(5)cpm=1-cpvcpd-1qcpd,where cpv=1870 J kg−1 K−1 is the isobaric heat capacity of dry air and cpd=1005 J kg−1 K−1 is the isobaric heat capacity of water vapor. Hence, cpm=1+0.86qcpd is a function of q, not a constant. The critical relative humidity rcontr can be calculated by G, Tcontr, and the ambient temperature T as follows:(6)rcontrT=GT-Tcontr+esatLTcontresatLT,where esatLT is the saturation pressure of water vapor over water and is given by:(7)esatLT=6.112exp⁡17.67T-273.15T-29.65,Following the Ponater et al. [20] study, we define a modified threshold rcrit∗ for contrail formation by combining rcontr with the threshold rcrit=0.7. Here rcrit is a critical threshold value that indicates a minimum relative humidity necessary to allow cloud formation within a grid box:(8)rcrit∗=rcontr·rcrit.because the PCC consists of SAC and the parameterization of natural cirrus coverage from temperature and humidity. In this paper, we adopt two different kinds of the parameterization of natural cirrus coverage to calculate and obtain two kinds of PCC fraction. One of the natural cirrus cloud coverage (bcirrus) from contrail parameterization in HadGEM2 is given by (Smith, 1990) the following:(9)bcirrusQn=0,ifQn≤-1121-Qn2,if0<Qn<11-121-Qn2,if0<Qn<11,if1≤Qn,where Qn is defined by(10)Qn=qiwcbcrit-1-r1-rcrit,where qiwc is the ice water content, bcrit=(1-rcrit)qsat(L), qsat(L) is the saturation specific humidity and r is the relative humidity over water. We employ rcrit∗ from (8) into (10) instead of rcrit and adopt the same bcrit: (11)Qn∗=qiwcbcrit-1-r1-rcrit∗.Meanwhile we also define and plunge Qn∗ into (9) instead of Qn. Therefore we obtain a potential cloud coverage for all high clouds (bcirrus(Qn∗)), which includes both nature cirrus and contrails. The first potential coverage for contrails called PCC1 could be calculated by(12)PCC1=bcirrusQn∗-bcirrusQn.Another the natural cirrus cloud coverage (bcirrus) from contrail parameterization in ECHAM4 is also given by (Sudqvist et al., 1989) the following:(13)bcirrusrcrit=0,ifrcrit≤r1-1-r-rcritrsat-rcrit,ifr<rcrit,where r is the relative humidity over ice in (13). Similarly, we can obtain bcirrus(rcrit∗) by replacing rcrit with rcrit∗ in (13). We use this parameterization of natural cirrus coverage to define a potential coverage for contrails called PCC2 as(14)PCC2=bcirrusrcrit∗-bcirrusrcrit.However, bcirrus calculated in this way is a three-dimensional field (longitude × latitude × level), so considering how to integrate a two-dimensional PCC field (longitude × latitude) over all 4 levels (from 200 hPa to 300 hPa), we adopt the maximum/random overlap assumption as described by Geleyn and Hollingsworth (1978).


## 3.1. Contrail Detection Method

Different from traditional CDA method and its variants only considering the spectral information of the image, our work focuses on the object-based method for the detecting contrails in each grid box of Himawari-8 images as shown in Figure2. The main process consists of the following two steps.


### 3.1.1. Multiresolution Segmentation

In our method, objects in Himawari-8 images are divided into different groups based on contrast, similarity, and shape of object, which are termed heterogeneity. There are several basic steps in the multiresolution segmentation procedure.(1)
The first step of segmentation procedure is to define a single pixel as the starting point.(2)
The single pixel fuses with other pixels as one group, if they are within the threshold of heterogeneity criteria.(3)
Pixel groups fuse with those neighbor groups.(4)
All of the pixels have been assigned to different objects until no further fusing is found.


### 3.1.2. Object-Based Classification

A nearest neighbor classification method is adopted to identify the Himawari-8 image into two classes: contrail and noncontrail object after multiresolution segmentation:(1)
It assigns several samples with contrail and noncontrail label, respectively. Objects should be verified with shape, size, and spectral characteristics.(2)
Contrail and noncontrail samples are trained from the segmented objects. Contrail and noncontrail classes are linked with the segmented image objects using the k-nearest neighbor classification method.(3)
Classification returns a fuzzy membership value for all the segmented image objects. A value of 0 means the object does not belong to contrail class, while a value of 1 means the object most likely belongs to contrail class.


## 3.1.1. Multiresolution Segmentation

In our method, objects in Himawari-8 images are divided into different groups based on contrast, similarity, and shape of object, which are termed heterogeneity. There are several basic steps in the multiresolution segmentation procedure.(1)
The first step of segmentation procedure is to define a single pixel as the starting point.(2)
The single pixel fuses with other pixels as one group, if they are within the threshold of heterogeneity criteria.(3)
Pixel groups fuse with those neighbor groups.(4)
All of the pixels have been assigned to different objects until no further fusing is found.


## 3.1.2. Object-Based Classification

A nearest neighbor classification method is adopted to identify the Himawari-8 image into two classes: contrail and noncontrail object after multiresolution segmentation:(1)
It assigns several samples with contrail and noncontrail label, respectively. Objects should be verified with shape, size, and spectral characteristics.(2)
Contrail and noncontrail samples are trained from the segmented objects. Contrail and noncontrail classes are linked with the segmented image objects using the k-nearest neighbor classification method.(3)
Classification returns a fuzzy membership value for all the segmented image objects. A value of 0 means the object does not belong to contrail class, while a value of 1 means the object most likely belongs to contrail class.


## 3.2. Artificial Visual Analysis

With the high temporal resolution stationary satellite images, we not only explore the occurrence of linear contrail but also study the development from linear contrail to contrail cirrus. In our evaluation, we investigate the satellite images of spring 2016 (March, May, and April) in our statistics analysis. After a large amount of observation analysis, we choose 30 min as the time scale of statistic, so as to avoid the overestimation of contrail, as well as fully considering the spreading of contrails. Thus, 4416 satellite images (24 per day and 92 days) are analyzed with both the proposed contrail detection method and artificial visual inspection. In the previous studies [1, 25–27], the reconstruction of contrail occurrence fraction adopts commercial flight track data to mask out regions of no air traffic. Their work has obvious limitations as follows:(1)
The polar satellite overpasses the observed region twice per day, which cause the severe insufficiency in observation sample.(2)
This kind of reconstruction did not consider the contrail spreading and persistence, with limited observation sample.Based on a very large amount of observation, some contrails could persist more than 5 h and spread 400–500 km at most, according to the annual mean wind speed at 200–300 hPa. Thus we define one new climatological index, termed contrail occurrence and persistence (COP). First, each of the grid boxes in satellite images is divided into 8 × 8 subgrid boxes; the COP is the ratio between the number of subgrids that contain detected contrails and the total number of subgrids, as illustrated in Figure4. Different from the traditional climatological index of contrail based on limited satellite overpass images, the COP index will count and aggregate detected contrails with 30 min time interval, which will reflect the occurrence and persistence of contrail.Figure 4
Illustration of8×8 subgrid box with detected linear contrails. The contrails occurrence fraction defined in this example is (8 + 8 + 4 + 3)/64 = 0.36.


## 3.3. Relative Humidity over Ice (RHI)

The contrails occurrence fraction presents a significant linear correlation versus the RHI on North America [25]. Therefore, for the comparison between the contrails occurrence fraction versus RHI on East Asia, we use the reanalysis data of temperature to calculate RHI by a formula as follows:(1)RHI=eesatiT≈qqsati,where q is the specific humidity, qsat(i) is the saturation specific humidity over ice, and esatiT is the saturation pressure of water vapor over ice and calculated by an empirical formula as follows:(2)esatiT=6.108exp⁡21.875T-276.16T-7.66,


## 3.4. Potential Contrail Coverage (PCC) Parameterization

Contrail formation depends on the surrounding temperature, moisture, and pressure fields that create an enough cold and moist ambient air for the mixing between the aircraft exhaust plume and its environment. Besides, we use SAC to diagnose the condition of ambient air with persistent contrails [16, 17]. The PCC is a method applied to determine the ability of the ambient air to allow contrail formation. An introduction necessary for PCC to estimate contrail coverage requires the fractional area of the atmosphere to support a persistent contrail formation [20]. Hence, in this paper, we adopt a similar method of PCC parameterization developed by Ponater et al. [20] for the ECHAM4 climate model. Relying on the thermodynamic theory, an empirical formula of the threshold values of temperature (Tcontr) and relative humidity (rcontr(T)) supporting contrail formation and persistence has been given according to Schumann [16]. In addition, contrail formation and persistence depend on some specific parameters and variables including the ambient temperature (T), pressure (p), specific humidity (q), the emission index of water vapor (EIH2o), the specific heat of fuel combustion (Q), and the propulsion efficiency of the aircraft (η). The thermodynamic theory predicts threshold values for temperature (Tcontr) and humidity (rcontrT). No contrails are possible if T>Tcontr and r<rcontrT. The formula of Tcontr (in °C) is given by the following:(3)Tcontr=-46.46+9.43ln⁡G-0.053+0.72ln⁡G-0.0532,where G is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute temperature versus water vapor partial pressure diagram and has the unit of Pa K−1, defined as follows:(4)G=EIH2o·cpmpεQ1-η,where EIH2o=1.25, ε=0.622 is the ratio of molecular masses of water and dry air, Q=4.3×107 J kg−1 K−1, η=0.3, and cpm is the isobaric heat capacity of moist air, which has the unit of J kg−1 K−1 and is given by the following:(5)cpm=1-cpvcpd-1qcpd,where cpv=1870 J kg−1 K−1 is the isobaric heat capacity of dry air and cpd=1005 J kg−1 K−1 is the isobaric heat capacity of water vapor. Hence, cpm=1+0.86qcpd is a function of q, not a constant. The critical relative humidity rcontr can be calculated by G, Tcontr, and the ambient temperature T as follows:(6)rcontrT=GT-Tcontr+esatLTcontresatLT,where esatLT is the saturation pressure of water vapor over water and is given by:(7)esatLT=6.112exp⁡17.67T-273.15T-29.65,Following the Ponater et al. [20] study, we define a modified threshold rcrit∗ for contrail formation by combining rcontr with the threshold rcrit=0.7. Here rcrit is a critical threshold value that indicates a minimum relative humidity necessary to allow cloud formation within a grid box:(8)rcrit∗=rcontr·rcrit.because the PCC consists of SAC and the parameterization of natural cirrus coverage from temperature and humidity. In this paper, we adopt two different kinds of the parameterization of natural cirrus coverage to calculate and obtain two kinds of PCC fraction. One of the natural cirrus cloud coverage (bcirrus) from contrail parameterization in HadGEM2 is given by (Smith, 1990) the following:(9)bcirrusQn=0,ifQn≤-1121-Qn2,if0<Qn<11-121-Qn2,if0<Qn<11,if1≤Qn,where Qn is defined by(10)Qn=qiwcbcrit-1-r1-rcrit,where qiwc is the ice water content, bcrit=(1-rcrit)qsat(L), qsat(L) is the saturation specific humidity and r is the relative humidity over water. We employ rcrit∗ from (8) into (10) instead of rcrit and adopt the same bcrit: (11)Qn∗=qiwcbcrit-1-r1-rcrit∗.Meanwhile we also define and plunge Qn∗ into (9) instead of Qn. Therefore we obtain a potential cloud coverage for all high clouds (bcirrus(Qn∗)), which includes both nature cirrus and contrails. The first potential coverage for contrails called PCC1 could be calculated by(12)PCC1=bcirrusQn∗-bcirrusQn.Another the natural cirrus cloud coverage (bcirrus) from contrail parameterization in ECHAM4 is also given by (Sudqvist et al., 1989) the following:(13)bcirrusrcrit=0,ifrcrit≤r1-1-r-rcritrsat-rcrit,ifr<rcrit,where r is the relative humidity over ice in (13). Similarly, we can obtain bcirrus(rcrit∗) by replacing rcrit with rcrit∗ in (13). We use this parameterization of natural cirrus coverage to define a potential coverage for contrails called PCC2 as(14)PCC2=bcirrusrcrit∗-bcirrusrcrit.However, bcirrus calculated in this way is a three-dimensional field (longitude × latitude × level), so considering how to integrate a two-dimensional PCC field (longitude × latitude) over all 4 levels (from 200 hPa to 300 hPa), we adopt the maximum/random overlap assumption as described by Geleyn and Hollingsworth (1978).


## 4. Result and Discussion

### 4.1. Contrails Detection

Figure5 shows the results of the detection using the proposed method in the different scenarios. The images were classified into two classes of contrail and noncontrail. The red area was classified as contrail and the blue area was classified as noncontrail. However, the contrails have still been broken into several parts and are disconnected with each other. As labeled in green cycle, there are some false alarms about contrail. Some misdetection happened and was labeled with the black cycle. The goal of contrail detection methods is to obtain a low false alarm rate (FAR) with sufficient detection efficiency (DE). Sensitivity of the algorithm depends on the thermal homogeneity and land/ocean contrast with contrail. Detection efficiency suffers from underestimation of contrail width and contrail cirrus overlap, while the algorithm cannot detect weak contrail structures that can still be seen by the human eye. Because of several uncertainties of algorithm identification, we also employ the artificial visual inspection to compute the COP fraction.Figure 5
The results of the object-based segmentation and classification under the three different scenarios: (a, b, c: left) is the original BTD image; (a, b, c: middle) is the result of segmentation; (a, b, c: right) is the final contrail classification result (green cycle: false alarms and black cycle: misdetection).
(a)
(b)
(c)


### 4.2. Plots of the COP Fraction versus the RHI and PCC

Figure6(a) presents the statistical COP fractions according to AI, while Figure 6(b) shows the statistical COP fractions derived from AVI. Although the general distribution of COP is similar to AVI, the AI field has obviously lower COP values than the AVI field because AI is hard to detect apotypic contrails like wide spread and fusion with natural cirrus. Due to the low spatial resolution, natural cirrus clouds with a linear appearance are also usually misclassified as contrails. The high rate of misdetections can be avoided by using higher spatial resolution satellite data. In our experiments, both the AI and AVI fields have broad maxima over east China’s Shandong, Anhui, and Jiangsu Province, minima over the sea of Japan and north of Korean Peninsula, another area of broad maxima around the Korean Strait and over the Yellow Sea, east of Japan’s Honshu Island. In addition, Figure 6(c) shows that the obvious difference between AI and AVI fields has broad large values over the region where AI and AVI fields have the similar large values. The reason for the difference should be that the linear or slightly bent shape of contrails is more easily identifiable by AVI. Although AI can identify more linear contrails over concentrated areas, the density of linear contrails identified by AVI is higher than AI over concentrated areas. It is reasonable that these obvious differences appear in the concentrated areas of linear contrails. However, we can reduce huge amount of statistical work through AI, and AVI has higher accuracy but shows powerless facing to massive satellite imagery processing, which is the reason that we only choose one season and a part of East Asia as the research time and area, respectively.Figure 6
COP fraction counted from Himawari-8 imagery by using AI (a) and AVI (b), and (c) presents the difference between (a) and (b) during spring (Mar, Apr, and May) 2016.
(a)
(b)
(c)As shown in Figure7, plots of the COP fractions versus the mean of RHI, PCC1, and PCC2 values calculated in 1°×1° grid box. Figures 7(a), 7(c), and 7(e) show the relationship between COP and variables for AI, and Figures 7(b), 7(d), and 7(f) show the same content for AVI. The Pearson correlation coefficients between COP and variables are presented in Table 3. As illustrated in Table 3 and Figure 7, there is a strong correlation between COP and RHI, corresponding to 0.72 for AVI, which is better than AI. And, it is obvious that the correlation coefficients for RHI are correspondingly greater than those for PCCs. In addition, PCC1’s linear correlations are clearly better than those for PCC2.Table 3
The Pearson correlation coefficient between COP and variables.
RHI
PCC1
PCC2
AI COP
0.59
∗
∗
0.52
∗
∗
0.45
∗
∗
AVI COP
0.72
∗
∗
0.68
∗
∗
0.39
∗
∗
Note. The double asterisks represent statistically significant values at 99% confidence level.Figure 7
Plot of COP fraction versus the mean of RHI, PCC1, and PCC2 values calculated in1°×1° grid box. (a), (c), and (e) represent AI. (b), (d), and (f) represent AVI.
(a)
(b)
(c)
(d)
(e)
(f)Contrails cause an alteration in the earth’s radiation budget, via their shortwave (SW) albedo effect and longwave (LW) greenhouse effect to impact climate. Contrail coverage and optical depth are essential aspects about contrail radiative forcing. Thus PCC1 (Smith 1990) reflect the actual reflection of contrail coverage better than PCC2 (Sudqvist et al., 1989) in model simulation, based on our observations and comparison.


## 4.1. Contrails Detection

Figure5 shows the results of the detection using the proposed method in the different scenarios. The images were classified into two classes of contrail and noncontrail. The red area was classified as contrail and the blue area was classified as noncontrail. However, the contrails have still been broken into several parts and are disconnected with each other. As labeled in green cycle, there are some false alarms about contrail. Some misdetection happened and was labeled with the black cycle. The goal of contrail detection methods is to obtain a low false alarm rate (FAR) with sufficient detection efficiency (DE). Sensitivity of the algorithm depends on the thermal homogeneity and land/ocean contrast with contrail. Detection efficiency suffers from underestimation of contrail width and contrail cirrus overlap, while the algorithm cannot detect weak contrail structures that can still be seen by the human eye. Because of several uncertainties of algorithm identification, we also employ the artificial visual inspection to compute the COP fraction.Figure 5
The results of the object-based segmentation and classification under the three different scenarios: (a, b, c: left) is the original BTD image; (a, b, c: middle) is the result of segmentation; (a, b, c: right) is the final contrail classification result (green cycle: false alarms and black cycle: misdetection).
(a)
(b)
(c)


## 4.2. Plots of the COP Fraction versus the RHI and PCC

Figure6(a) presents the statistical COP fractions according to AI, while Figure 6(b) shows the statistical COP fractions derived from AVI. Although the general distribution of COP is similar to AVI, the AI field has obviously lower COP values than the AVI field because AI is hard to detect apotypic contrails like wide spread and fusion with natural cirrus. Due to the low spatial resolution, natural cirrus clouds with a linear appearance are also usually misclassified as contrails. The high rate of misdetections can be avoided by using higher spatial resolution satellite data. In our experiments, both the AI and AVI fields have broad maxima over east China’s Shandong, Anhui, and Jiangsu Province, minima over the sea of Japan and north of Korean Peninsula, another area of broad maxima around the Korean Strait and over the Yellow Sea, east of Japan’s Honshu Island. In addition, Figure 6(c) shows that the obvious difference between AI and AVI fields has broad large values over the region where AI and AVI fields have the similar large values. The reason for the difference should be that the linear or slightly bent shape of contrails is more easily identifiable by AVI. Although AI can identify more linear contrails over concentrated areas, the density of linear contrails identified by AVI is higher than AI over concentrated areas. It is reasonable that these obvious differences appear in the concentrated areas of linear contrails. However, we can reduce huge amount of statistical work through AI, and AVI has higher accuracy but shows powerless facing to massive satellite imagery processing, which is the reason that we only choose one season and a part of East Asia as the research time and area, respectively.Figure 6
COP fraction counted from Himawari-8 imagery by using AI (a) and AVI (b), and (c) presents the difference between (a) and (b) during spring (Mar, Apr, and May) 2016.
(a)
(b)
(c)As shown in Figure7, plots of the COP fractions versus the mean of RHI, PCC1, and PCC2 values calculated in 1°×1° grid box. Figures 7(a), 7(c), and 7(e) show the relationship between COP and variables for AI, and Figures 7(b), 7(d), and 7(f) show the same content for AVI. The Pearson correlation coefficients between COP and variables are presented in Table 3. As illustrated in Table 3 and Figure 7, there is a strong correlation between COP and RHI, corresponding to 0.72 for AVI, which is better than AI. And, it is obvious that the correlation coefficients for RHI are correspondingly greater than those for PCCs. In addition, PCC1’s linear correlations are clearly better than those for PCC2.Table 3
The Pearson correlation coefficient between COP and variables.
RHI
PCC1
PCC2
AI COP
0.59
∗
∗
0.52
∗
∗
0.45
∗
∗
AVI COP
0.72
∗
∗
0.68
∗
∗
0.39
∗
∗
Note. The double asterisks represent statistically significant values at 99% confidence level.Figure 7
Plot of COP fraction versus the mean of RHI, PCC1, and PCC2 values calculated in1°×1° grid box. (a), (c), and (e) represent AI. (b), (d), and (f) represent AVI.
(a)
(b)
(c)
(d)
(e)
(f)Contrails cause an alteration in the earth’s radiation budget, via their shortwave (SW) albedo effect and longwave (LW) greenhouse effect to impact climate. Contrail coverage and optical depth are essential aspects about contrail radiative forcing. Thus PCC1 (Smith 1990) reflect the actual reflection of contrail coverage better than PCC2 (Sudqvist et al., 1989) in model simulation, based on our observations and comparison.


## 5. Conclusion

Contrails are important for understanding how the humans contribute to climate change by one artificial element into the atmosphere. Identifying contrails is the first step in contrail relative study. In this paper, object-based contrail detecting method is proposed to reflect the process of contrail occurrence and persistence. According to the algorithm identification (AI) and artificial visual inspection (AVI), COP measured from Himawari-8 stationary satellite imagery is related to upper tropospheric relative humidity over ice (RHI) computed with the ECMWF reanalysis data by simple linear regression. Two kinds of potential contrail coverage (PCC1 and PCC2), from contrail parameterization in ECHAM4 and HadGEM2, are compared with distribution of COP. Through the linear correlation between COP and PCCs fractions, we found the PCC1 has better correspondence with COP, which means that the PCC1 in model simulation reflects the distribution of actual contrails better than PCC2.



---

*Source: 8707234-2017-10-01.xml*