Northern Hemisphere Contrail Properties Derived from Terra and Aqua MODIS Data for 2006 and 2012

Linear contrail coverage, optical property, and radiative forcing data over the Northern Hemisphere (NH) are derived from a year (2012) of Terra and Aqua Moderate-resolution Imaging Spectroradiometer (MODIS) imagery, and are compared with previously published 2006 results (Duda et al., 2013; Bedka et al., 2013; Spangenberg et al., 2013) using a 10 consistent retrieval methodology. Differences in the observed Terra-minus-Aqua screened contrail coverage and patterns in the 2012 annual-mean air traffic estimated with respect to satellite overpass time suggest that most contrails detected by the contrail detection algorithm (CDA) form approximately 2 h before overpass time. The 2012 screened NH contrail coverage (Mask B) shows a relative 3% increase (from 0.136% to 0.140%) compared to 2006 data for Terra and increased by almost 7% (0.134% to 0.143%) for Aqua. A new post-processing algorithm added to the contrail mask processing estimated that the 15 total contrail cirrus coverage visible in the MODIS imagery may be three to four times larger than the linear contrail coverage detected by the CDA. This estimate is similar in magnitude to the spreading factor estimated by Minnis et al. (2013). Contrail property retrievals of the 2012 data indicate that both contrail optical depth and contrail effective diameter decreased approximately 10% between 2006 and 2012. The decreases may be attributed to better background cloudiness characterization, changes in the waypoint screening, or changes in contrail temperature. The total mean contrail radiative 20 forcing (TCRF) for all 2012 Terra observations were -6.3, 14.3, and 8.0 mW m for the shortwave (SWCRF), longwave (LWCRF), and net forcings, respectively. These values are approximately 20% less than the corresponding 2006 Terra estimates. The decline in TCRF results from the decrease in normalized CRF, partially offset by the 3% increase in overall contrail coverage in 2012. The TCRFs for 2012 Aqua are similar, -6.4, 15.5, and 9.0 mW m for shortwave, longwave, and net radiative forcing. The strong correlation between the relative changes in both total SWCRF and LWCRF between 2006 25 and 2012 and the corresponding relative changes in screened contrail coverage over each air traffic region suggests that regional changes in TCRF from year to year are dominated by interannual changes in contrail coverage over each area.

1 Introduction Duda et al. (2013) modified the Mannstein et al. (1999) contrail detection algorithm (CDA) to use Moderate Resolution Imaging Spectroradiometer (MODIS) brightness temperature difference (BTD) data.The CDA was part of the development of a consistent global observational dataset of linear contrail properties including coverage, optical properties (Bedka et al., 2013), and radiative forcing (Spangenberg et al., 2013) from a single satellite imager analyzed with the same methodology.Such a dataset is useful for validating and improving the representation of contrail cirrus within climate models (Burkhardt et al., 2010).We evaluate the accuracy of the CDA from Duda et al. by using a visual analysis tool (Minnis et al., 2005) to allow human analysts to subjectively determine linear contrail coverage in satellite imagery and by comparing the analyst coverage with the coverage determined by the CDA.The CDA and the visual analysis procedure are briefly described in section 2, and the results of the visual analysis, including the determination of a corrected coverage based on an evaluation of the CDA's overall detection efficiency, are presented in section 3. Duda et al. (2013) developed a modified version of the Mannstein et al. (1999) technique to determine the linear contrail coverage over the Northern Hemisphere (NH) during 2006 from Aqua and Terra MODIS radiance imagery.The original CDA was designed for AVHRR (Advanced Very High Resolution Radiometer) data, and when applied to MODIS radiances, the number of false detections increased significantly even when the sensitivity of the algorithm was adjusted to account for the greater sensitivity of the MODIS to cloud structures in the IR imagery.To reduce the number of false detections, other thermal infrared channels (6.7, 8.6, and 13.3 µm) available on MODIS were used to eliminate the detection of surface and lower cloud features as contrails, and global aircraft emissions waypoint (AEDT) data provided by the FAA (Wilkerson et al., 2010) allowed comparison of detected contrails with commercial aircraft flight tracks.

Contrail detection algorithm
Like the Mannstein et al. algorithm, the Duda et al.CDA uses only brightness temperature (BT) data from the IR channels of MODIS and can be applied to both day and night scenes.A scene-invariant threshold is used to detect cloud edges produced by contrails, and a set of 4 binary masks determines if the detected linear features are truly contrails.This modified CDA follows the same data flow as the original Mannstein et al.CDA.The IR brightness temperature measurements are first inverted (multiplied by -1) so they can be treated with the same procedures as the BTD measurements (for example, T12 is the inverted BT for the 12.0 μm measurements).To normalize the brightness temperature (and BTD) data, the BT data are smoothed with a rotationally symmetric Gaussian 5×5 pixel lowpass kernel (resulting in, for example, T12 $$$$$ ), and the local standard deviation of the data are computed over the region of the kernel by filtering the square of the differences between the original and smoothed images again with the same kernel: (1) The difference between the original and the smoothed images was then normalized with the local standard deviation to get the normalized image: (2) where 0.1 K is added to the local standard deviation in (2) to limit the sensitivity of the normalization over very homogeneous regions.In addition to T12 and the 11 -12 μm BTD (hereafter called BTD1) data used in Mannstein et al. (1999), the modified CDA also uses 6.8 μm water vapor channel BT (T6.8), the BTD between 8.6 and 12.0 μm (BTD2), and the BTD between 8.6 and 13.3 μm (BTD3).A fourth brightness temperature difference was computed from a sum of the BTDs (BTD4 = BTD1 + BTD2 -T6.8).
The sum of the normalized images N = N12+NBTD1+NBTD3 was then convolved with a line filter of 19×19 pixels in 16 different directions, and the individual connected regions resulting from the filtering are considered as possible contrail objects.These objects are then checked against 4 binary masks to check for contrails.The first three masks are similar to the original CDA: The threshold in condition AA selects pixels having a higher than average normalized brightness, which is typical of thin cirrus in T12 or BTD imagery.The low BTD1 threshold of the original CDA was retained in condition BB to allow detection of possible contrails over lower opaque clouds, but an upper threshold of 4.5 K was added to remove some liquid water clouds that have an exceptionally large BTD1.Contrails are usually not optically thick enough to produce BTD1 as  regions (e.g., United States, Europe) while removing false detections in low traffic regions (Duda et al., 2013).The flight track screening is especially useful in the tropics where thin cirrus streaks from convection are often misinterpreted by the CDA as linear contrails.It does not, however, account for military flights which are not included in the flight track database.

Visual analysis program
To test the accuracy of the CDA, a team of 5 human analysts visually evaluated a group of 40 satellite images (5minute MODIS granules), and the results were compared with the output from the CDA.The sample size was limited due to the intensive labor requirements of the visual analysis.To create a representative set of imagery, the satellite granules were randomly selected from the 2006 Aqua and Terra MODIS data.The selected granules were split evenly (20 each) between Aqua and Terra, and between daytime and nighttime scenes.In addition, the granules were selected to encompass the range of contrail numbers detected per image (from near zero to several thousand), so that the selected granules would accurately represent the overall distribution of contrail coverage (number of contrail pixels per granule) determined from Mask B from the 2006 MODIS data.Eight granules each were randomly selected from five contrail number intervals (Table 2) so that the selected granules would represent the overall contrail density distribution.(The selected contrail number intervals comprise 99.4 percent of all granules and 96.5 percent of detected contrail pixels.) Figure 1 presents the locations of the selected Aqua and Terra granules, and Table 3 shows the time and date of each granule selected for visual analysis.Most granules occur over high air traffic regions (CONUS, Europe, eastern Asia), and correspond to the regions where most contrails were detected in the satellite imagery.
To determine the overall accuracy of the contrail masks, the subjective error analysis was performed using an updated version of the interactive program developed by Minnis et al. (2005).As in Minnis et al., the visual analysis program displays both T11 and BTD1 images with removable overlays of the CDA-detected contrails.The analyst can evaluate the images using a range of contrast adjustments, and add or delete pixels based on the analyst's judgment.The analyst can either add or subtract contrail pixels by selecting rectangular boxes within a magnified (zoom) window.Additionally, visible imagery (0.65 μm) and multi-channel RGB imagery were also used (when available) to help assess whether the detected contrails were accurate or not.(During the daytime, the RGB image was a combination of visible, 3.7 -11 μm and 11 μm imagery, while it was a combination of 3.7 μm, 11 μm, and 3.7 -11 μm imagery at night.)Finally, commercial flight track information was provided in the visual analysis program through imagery that showed flight paths within the satellite field of view up to 4 h before overpass time.The five analysts interactively evaluated Mask B images to estimate the number of false detections (deleted pixels), missed contrails (added pixels), and positive contrail detections (retained pixels).A composite mask for each image was determined from a majority consensus of the five analysts: a contrail pixel was retained or deleted, or new pixels added, if at least three of the analysts agreed with such an assessment.Because of the subjective nature of the analysis, the consensus choice of the analysts is expected to be superior to any single visual analysis.
Figure 2 shows an example of the process for a Terra MODIS scene.The BTD image (Figure 2a) reveals many contrails in the image along with some apparent contrail cirrus.Mask A (shown as red lines in Figure 2b) coincides with some of the contrails.Mask B (yellow) adds to Mask A, lengthening and widening the results to give more coverage.Mask C (green) appears to detect nearly every contrail and more coverage.The same masks also identify the lower right linear features as contrails.The analysts added (green in Figure 2c) many pixels to Mask A (red), and deleted (blue) all of the contrail pixels in the lower right.The added and retained pixels then serve as the best reference mask for evaluating the other two masks.For this case, both B and C correctly identified many more contrail pixels than A, but also detected more false contrail pixels.The cloud streets in the lower right corner of the images are liquid water clouds with a BTD1 signature similar to the contrails in the image center.Mask C detects nearly all of these clouds as contrails, due to their relative orientation and proximity to the flight tracks (not shown) in the region.

Visual analysis of contrail coverage
For all analyzed images, Mask C had the most false alarms, sometimes identifying cloud features perpendicular to the flight track as contrails, while masks A and B generally identified contrails that were parallel to the flight track orientation.Mask B, with the exception of a few cases of very thin contrails, also identified most of the linear contrails correctly identified by Mask C. Mask B correctly detected the main body, but missed the endpoints, of the very thin contrails.Mask A detected some of the linear contrails, but missed many spreading contrails as well as the contrail endpoints as in Mask B. Additionally, Mask A frequently missed linear contrails when the BTD was small.The bias ratio ((retained+deleted) pixels/(retained+added) pixels) for masks A, B, and C were 0.67, 1.32 and 4.55 respectively.Based on the relatively equal numbers of deleted and added pixels from the visual analysis suggests that Mask B is the most appropriate mask for measuring linear contrail coverage.As shown below, the relative number of false detections may not be independent of location and may be more common in Southeast Asia, but a larger sample size of analyzed granules would be necessary to confirm any geographic variation in false positive detections.
While it is assumed that the analysts' subjective assessments represent the "ground truth" set for the CDA selection, subjectively defining a linear contrail is not always straightforward.This difficulty arises from several factors including the presence of cirrus and older contrails that have spread significantly.When it was not possible to decide if the pixels were actually contrails or cirrus, the CDA result was not changed.The most common source of misdetections was cirrus streamers from mid-latitude storms and tropical convection.Underestimates of contrail coverage were most common within clusters of old, overlapping and decaying contrails where the appearance of the contrails has changed considerably from a typical newly formed contrail.Assessment of these regions by the visual analysts was especially difficult.No explicit attempt was made to add all contrail pixels from spreading contrails due to the subjective nature of the cloudiness.Multiple analysts were used to help minimize some of the uncertainty in contrail cover, but this evolution in contrail appearance is common in contrail outbreaks (Minnis et al., 1998;Carleton et al., 2008).Such outbreaks may be responsible for much of the excess cloudiness due to contrails, but are difficult to quantify with automated methods.The ambiguity resulting from older and overlapping contrails in such outbreaks should be considered when defining linear contrail extent (Minnis et al., 2013).
Thus, it is likely that the truth set from the analysts is an underestimate and the corresponding CDA results will also underestimate the true linear contrail and contrail cirrus coverage.The FAR is defined as the ratio of the pixels falsely classified as contrails to the total number of analyzed pixels, and the DEF is the ratio of the contrail pixels correctly identified by the CDA to the total number of true contrail pixels (identified from the visual analysis) (Meyer et al. (2002)).Figure 3 presents the false alarm rate of contrail detection computed from the visual analysis for the Aqua MODIS granules.The results are roughly comparable to Meyer et al. (2002).The similarity is likely due to most of the Aqua granules being centered over the high air traffic regions over the mid-latitudes (eastern CONUS, north Atlantic corridor, and Europe).No significant relationship between the false alarm rate and the 12-μm BT standard deviation was found for the Terra MODIS granules, as the range in 12-μm BT standard deviation was narrow, and several of the Terra visual analysis granules are located in southeast Asia where the false alarm rate is on average 0.06% higher (0.14% versus 0.08%) than the Europe/Northern Atlantic/CONUS granules.The Terra results suggest that FAR may vary across the Northern Hemisphere, but as most air traffic occurs over Europe/Northern Atlantic/CONUS, the Aqua-based FAR are considered to be the most representative of the FAR -SDT 12 $$$$$$$$ relationship over the NH.Following Meyer et al. (2002), the false alarm corrected contrail frequency ( )* ++ ) is computed from the uncorrected contrail coverage ( )* + ) by: three-parameter contrail correction methodology to the MODIS data for the range of air traffic amounts.We assume that the relation between SDT 12 $$$$$$$$ and the false alarm-corrected contrail coverage is constant in each air traffic intensity region, and SDT 12 $$$$$$$$ is not correlated with air traffic density.Regions with an annual mean SDT 12 $$$$$$$$ greater than 1.2 K (parts of the Himalayas) were not analyzed due to the large background inhomogeneity and the small number of contrails detected there.

Determination of background heterogeneity-corrected coverage
Contrail coverage in each 1×1 degree grid box was segregated into five air traffic density bins (Table 4).The air traffic density for each grid box was estimated from the number of flight segments above 7.15 km tabulated during the 4-h period before satellite overpass time on 15 October 2006.The bin with the lowest air traffic (low) accounts for 80 percent of all grid boxes (no air traffic is found in roughly half of all grid boxes -because so many grid boxes had no air traffic we use only grid boxes with at least 10 flight segments per grid box (low4)), the "low-mid" bin accounts for the 80 through 90 th percentile of air traffic, the "mid" bin covers the 90 to 95 th percentile of the NH air traffic distribution.The "mid high" bin accounts for the grid boxes with the 95 to 99 percentile of air traffic and the final bin (high) covers the top percentile in air traffic density.Figure 4 We estimated the contrail detection efficiency visually from the 40 MODIS granules.Figure 5 shows the contrail detection efficiency computed from the visual analysis granules.Unlike Meyer et al. (2002), we compute DEF as a function of SDT 12 $$$$$$$$ : DEF = 0.785 -0.155 SDT 12 $$$$$$$$ .The mean and standard deviation of DEF for all 40 visually analyzed granules is 0.68 ± 0.23 (compared to 0.4± 0.2 in Meyer et al. (2002)).
From the false alarm corrected contrail frequency, the detection efficiency, and the relationships between contrail coverage and SDT 12 $$$$$$$$ for each air traffic density range, the corrected contrail coverage was computed for both Aqua and Terra during 2006 from equations ( 3), (4), and (5). Figure 6 shows the NH annual-mean corrected contrail coverage for both satellites during 2006.The corrected coverage increases 25% for Terra (from 0.136 to 0.170%), while the Aqua coverage increased from 0.134 to 0.169%), a relative increase of 26%.The corrected Terra coverage shows that the maxima in linear contrail coverage are now over CONUS and Europe, although the corrected Aqua coverage still has a maximum over the North Atlantic due to the decrease in detectability of linear contrails over the continents during the afternoon after several hours of air traffic have potentially saturated the skies with multiple contrails.

Analysis of retrieved contrail properties
Following the completion of the contrail coverage analysis with the CDA, optical properties of the detected contrails were computed (Bedka et al., 2013).The retrieved contrail optical properties (optical depth (OD), effective particle size (De)) did not change significantly between the retained, added, and deleted contrails in the 40 analyzed granules.This result is expected as the binary masks are designed to detect primarily thin ice clouds with a radiative signature similar to linear contrails.
A geographical analysis of the contrail optical properties showed that the retrieved optical depths and particle sizes in the contrails were generally consistent throughout the Northern Hemisphere except for some regions in the Arctic where the number of acceptable retrievals was small, and over the Tibetan plateau, where retrieved contrail optical depths appeared to be anomalously large, especially during April 2006.Fifty-nine MODIS granules (29 Aqua, 30 Terra) from April 2006 passing over Tibet were visually examined by one of the analysts to determine nine (4 Aqua, 5 Terra) samples of images where the contrail detection was particularly poor, and nine samples where the contrail detection was relatively accurate.
The sample of MODIS granules showed that the retrieved properties in the poor detection cases (mean OD = 0.52, mean De = 48.7 μm) were significantly different from the properties in the better detection cases (mean OD = 0.34, mean De = 38.4μm).Filtering out pixels where the standard deviation in the 12-μm BT was much larger than normal can eliminate some contrail misdetections, but the regional average properties over Tibet were still a bit larger than the hemispheric average for April 2006 (mean OD = 0.20, mean De = 36.7 μm).
In addition to the binary masks, other quality control methods including de-striping and remapping of the MODIS imagery to a standard projection, and checking flight track data to confirm contrail formation were used to reduce the number of false positive detections in the CDA.In particular, global aircraft emissions waypoint (AEDT) data provided by the FAA(Wilkerson et al., 2010) allowed comparison of detected contrails with commercial aircraft flight tracks.A pixellevel product based on the advected flight tracks defined by the waypoint data and the three-dimensional wind component profiles from the Goddard Earth Observing System (GEOS-4) reanalyses(Bloom et al., 2005) was developed to assign a confidence of contrail detection (CCD) for the contrail mask pixels.The product computes a flight track "strength" value for each pixel in a MODIS granule determined from the relative age of the most recent advected flight tracks passing over the pixel.Flight tracks as old as 4 hours before image time are considered because the young contrails developing from these tracks are most likely to be detected in the MODIS imagery(Kärcher et al., 2010).To compute the CCD parameter, a subjective error analysis was performed using an updated version of the interactive program developed byMinnis et al. (2005) (described in more detail in section 3).The analysis was applied to sets of 22 (Terra and Aqua) MODIS images over the North Atlantic flight corridor as well as a few Terra MODIS images over North America.The corridor was selected as it has large amounts of air traffic with a relatively simple zonal flight track orientation to allow for easier visual detection of linear contrails.Four different analysts interactively evaluated the Mask B images to estimate the number of false detections (deleted pixels), missed contrails (added pixels), and positive contrail detections (retained pixels) for the contrail mask.A composite mask for each image was determined by consensus of the four analysts.The visual analyses served as "ground truth" within a logistic regression model to compute the likelihood that a CDA contrail pixel was actually a contrail.The logistic regression model uses several variables to determine the CCD of each contrail pixel in the composite mask.The variables in the model include the number of advected flight tracks passing over the region, the pixel's flight track strength value, and the bearing of the flight track relative to the orientation of the detected contrail, which are used to compute the CCD for each contrail pixel.The CCD is used to filter out pixels that do not correspond to any known flights, thus reducing the false detections in all three masks.The flight track screening retains nearly all of the contrails detected in the high traffic Mannstein et al. (1999) andMeyer et al. (2002) both noted that contrail detection by the CDA is affected by the background heterogeneity of the thermal IR imagery.The detection efficiency is highest over homogeneous backgrounds (cloudless ocean, flat land) and decreases over heterogeneous backgrounds (mountains, broken cloudiness).Mannstein et al.     developed a one-parameter correction for background heterogeneity from a linear regression of contrail coverage with the standard deviation of the 12-μm BT, while Meyer et al. developed a more advanced three-parameter correction based on the false alarm rate (FAR), the detection efficiency (DEF) of the CDA, and the heterogeneity of the background in the thermal IR imagery, quantified here as the local standard deviation of the 12-μm BT (SDT 12 $$$$$$$$ ).
While the Meyer et al. study used AVHRR data over central Europe, a region with high air-traffic densities, the Duda et al. study surveyed the Northern Hemisphere with a wide range of air traffic densities.We applied the Meyer et al.
shows the contrail coverage as a function of the standard deviation of the 12-μm BT for each of the air traffic density bins.Comparing from the low air traffic to the high air traffic plot, both the slope and the correlation coefficient of the regression line increase as air traffic density increases.The increase in slope indicates that the heterogeneity of the background becomes more important as air traffic density increases and the largest corrections are needed in the regions with the highest air traffic density.The regressions for each air traffic region are used to compute the homogenized contrail coverage : ct,min A: the slope and In is the intercept of the air traffic-sampled regressions in Figure4.As stated inMeyer et al. (2002), the homogenized contrail coverage  ct,min represents a lower limit for the true area covered by linear contrails because the DEF of the CDA is low due to the tuning of the algorithm to a low FAR.The resulting corrected contrail coverage  )

Figure 1 .Figure
Figure 1.Location of Aqua (a) and Terra (b) MODIS granules selected for visual analysis.

Figure 3 .
Figure 3. False alarm rate of contrail detection computed from visual analysis.Regression line computed from averaged results represented by the four plus signs in the plot.

Figure 4 .
Figure 4. Contrail coverage determined from 2006 Aqua and Terra MODIS imagery versus the standard deviation of the 12 µm BT for five air traffic densities.

Figure 5 .
Figure 5. Contrail detection efficiency computed from visual analysis.Regression line computed from averaged results represented by the four plus signs in the plot.

Figure 6 :
Figure 6: The 2006 annual-mean NH contrail coverage corrected for background inhomogeneity in Terra and Aqua MODIS imagery.