Ash rich clouds, produced by explosive volcanic eruptions, are a major hazard to aviation. Unfortunately, explosive volcanic eruptions are not always detected in a timely manner in satellite data. The large optical depth of emergent volcanic clouds greatly limits the effectiveness of multispectral infrared‐based techniques for distinguishing between volcanic and nonvolcanic clouds. Shortwave radiation‐based techniques require sufficient sunlight and large amounts of volcanic ash, relative to hydrometeors, to be effective. Given these fundamental limitations, a new automated technique for detecting emergent clouds, produced by explosive volcanic eruptions, has been developed. The Cloud Growth Anomaly (CGA) technique utilizes geostationary satellite data to identify cloud objects, near volcanoes, that are growing rapidly in the vertical relative to clouds that formed through meteorological processes. Explosive volcanic events are shown to frequently be a source of rapidly developing clouds that, at a minimum, reach the upper troposphere. As such, the CGA algorithm is effective at determining when a recently formed cloud is possibly the result of an explosive eruptive event. While the CGA technique can be applied to any geostationary satellite sensor, it is most effective when applied to latest generation of meteorological satellites, which provide more frequent images with improved spatial resolution. Using a large collection of geographically diverse explosive eruptions, and several geostationary satellites, the CGA technique is described and demonstrated. A CGA‐based eruption alerting tool, which is designed to improve the timeliness of volcanic ash advisories, is also described.

1 Introduction Ash‐rich clouds produced by explosive volcanic eruptions are a major hazard to aviation. Flying in volcanic ash‐contaminated air can cause severe damage to jet engines and airframes and, in the worst case, cause in‐flight engine failure (Casadevall, 1994; Clarkson et al., 2016; Guffanti et al., 2010; Miller & Casadevall, 2000). As documented by Guffanti et al. (2010), the most severe (damaging) encounters with volcanic ash have generally occurred within 1,000 km of the source volcano in the early stages of volcanic constituent transport and dispersal in the atmosphere. In addition, volcanic clouds produced by explosive eruptions can reach jet aircraft cruising altitudes in as little as 5 min and 10 or more eruptions per year will generate volcanic clouds with tops at or above jet cruising altitudes (International Civil Aviation Organization, 2007). Thus, timely and accurate detection of explosive volcanic activity is critical for aviation safety. Unfortunately, explosive volcanic eruptions are not always detected in a timely manner, especially at unmonitored volcanoes. Approximately 90% of the world's volcanoes are not regularly monitored for activity (Ph. Bally Ed, 2012). For instance, when Nabro volcano in Eritrea erupted for the first time in recorded history on 12 June 2011 (Global Volcanism Program, 2011), the eruption went undetected by authorities for several hours because it occurred at a remote, unmonitored, volcano with no past history (Wiart & Oppenheimer, 2005) and the resulting cloud was difficult to distinguish from meteorological clouds in satellite imagery. Volcanic Ash Advisory Centers (VAACs) are responsible for issuing volcanic ash advisories for aviation. While VAAC forecasters utilize an array of tools to detect new ash emissions and track existing ash clouds, meteorological satellites are the primary tools utilized by VAACs (e.g., Prata, 2009; Tupper et al., 2004). Geostationary satellites, in particular, are most often used to detect explosive volcanic eruptions. Volcanic clouds formed as a result of explosive volcanism evolve rapidly in time, so the frequent images provided by geostationary satellites are essential for capturing the early stages of an eruption. Geostationary satellites such as Himawari‐8 (launched in 2014; Bessho et al., 2016) and the next generation Geostationary Operational Environmental Satellite (GOES‐R; Schmit et al., 2017) provide much more frequent imagery, with 4 times the spatial resolution, compared to the previous generation of satellites. The imaging radiometers on the Himawari‐8 and GOES‐R satellites also have more spectral channels than the previous generation. Despite the improved spectral, spatial, and temporal capabilities of geostationary satellites such as Himawari‐8, timely volcanic eruption detection remains a challenge. 2006 2015b 2001 2000 2004 1997 2006 2007 Explosive volcanic eruptions can be difficult to manually distinguish from meteorological clouds in convectively active environments, especially early in the evolution. From a practical standpoint, it is difficult to carefully manually analyze all satellite images over all volcanoes in near real time. The data volume issue is becoming more pronounced as next generation geostationary satellites, such as Himawari‐8, produce approximately 100 times more data each day than the previous generation of geostationary satellites. The large optical depth of emergent volcanic clouds greatly limits the effectiveness of multispectral infrared‐based techniques for detecting volcanic ash from satellites (Pavolonis et al.,; Prata et al.,; Simpson et al.,; Tupper et al.,). When sufficient sunlight is present and volcanic ash is not encased in ice, ultraviolet (Seftor et al.,) and visible (Pavolonis et al.,) radiation‐based techniques are sometimes useful for detecting ash in clouds that are opaque to infrared radiation, but such techniques cannot be relied on for all eruptions. Thus, VAAC forecasters often utilize manual analysis of time sequences of satellite imagery to extract features that exhibit volcanic cloud like attributes such as sudden cloud development, over a volcano, which is out of sync with the surrounding meteorology (Tupper et al.,). While the human expert analysis techniques are effective, there are two important limitations. Given the limitations above, the goal of this manuscript is to introduce a new automated technique that utilizes infrared‐based cloud vertical growth rates (e.g., trends in metrics related to cloud top temperature) to detect explosive volcanic eruptions. The cloud growth‐based methodology, referred to as the Cloud Growth Anomaly (CGA) technique, can be applied to any geostationary satellite (day and night) but is particularly valuable when applied to next generation satellites such as Himawari‐8 and GOES‐R. The CGA technique utilizes infrared measurements to identify cloud objects and compute cloud vertical growth rates from two successive images. As will be shown, growth rates, associated with a large collection of volcanic eruptions, are nearly always statistically distinguishable from meteorological cloud growth rates. In the remainder of this paper, the cloud object identification and tracking technique will be described. In addition, the vertical growth rates associated with 79 explosive volcanic events, from around the world, are compared to a very large collection of meteorological clouds imaged by several geostationary satellites with widely varying spatial, spectral, and temporal capabilities. Near‐real‐time applications of the CGA technique, such as volcanic eruption alerting and mass eruption rate determination, will also be discussed.

2 Cloud Object Identification A recursive procedure was developed to extract cloud objects from infrared satellite imagery in a manner that is reasonably consistent with human expert identification of localized cloud elements. As in Cintineo et al. (2013, 2014), cloud objects are extracted from the 11‐μm top of troposphere cloud emissivity parameter (ε tot ) first introduced by Pavolonis (2010). The ε tot is simply the 11‐μm effective cloud emissivity a cloud would have if the radiative center were located at the top of the troposphere. Every meteorological satellite in geostationary orbit (past, present, and future) has an imaging radiometer with a spectral channel centered near 11 μm, so the cloud object identification and tracking methodology introduced in this manuscript can be applied to any geostationary satellite. Unlike the 11‐μm brightness temperature (BT), ε tot is constrained to range between 0.0 and 1.0, where confidence that a cloud is present greatly increases as ε tot increases. Examples of the ε tot parameter and the corresponding false color ash/dust image (Lensky & Rosenfeld, 2008; Pavolonis et al., 2015a) and the 11‐μm brightness temperature image are shown in Figure 1. The first column in Figure 1 is derived from a Spinning Enhanced Visible and InfraRed Imager (SEVIRI) image taken at 20:30 UTC on 12 June 2011, and the second column was derived from the next SEVIRI image at 20:45 UTC. Note in Figure 1 (bottom row) how regions of localized maxima in the ε tot image can be readily identified and correspond to localized minima in the 11‐μm BT image. The goal of the cloud object identification procedure is to automate the extraction of cloud object features that correspond to regions of localized maxima (minima) in the ε tot (11‐μm BT) image, consistent with the visual appearance of an opaque, near‐source, volcanic cloud in infrared satellite imagery (e.g., Schneider & Rose, 1994; Tupper & Kinoshita, 2003). Figure 1 Open in figure viewer PowerPoint (first row) SEVIRI false color imagery, (second row) 11‐μm brightness temperature (BT), and (third row) 11‐μm top of troposphere emissivity (ε tot ) are shown for two consecutive image times. The first column is from the 20:30 UTC Meteosat‐9 (MSG‐2) SEVIRI scan on 12 June 2011. The second column is from the 20:45 UTC Meteosat‐9 SEVIRI scan on 12 June 2011. The perimeter of automatically defined cloud objects is overlaid on each image in black. Using the same methodology as Wielicki and Welch (1986) and Pavolonis et al. (2015b), an initial set of parent cloud objects is constructed using all pixels where ε tot is greater than or equal to 0.05. Each parent cloud object is recursively decomposed into smaller objects by incrementally increasing the ε tot threshold for object membership (objects are constructed from pixels with ε tot ≥ the threshold value) until the use of larger ε tot threshold values no longer results in more than one object. The ε tot cloud object membership threshold values utilized in sequential order are 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.75, 0.80, 0.85, and 0.90. As Figure 1 shows, the cloud object identification algorithm is able to capture regions of localized maxima in ε tot in a manner that is generally consistent with human pattern recognition. In particular, the boundaries of the emergent Nabro cloud from the 12 June 2011 eruption (Global Volcanism Program, 2011) are consistent with the infrared satellite imagery. As will be shown using additional examples, the recursive cloud object identification methodology is able to routinely determine the boundaries of emergent volcanic clouds in a manner that is consistent with human interpretation of satellite imagery.

3 Cloud Object Tracking In order to determine if an existing cloud object has undergone vertical growth or if a new object has emerged from the background, cloud objects must be tracked from one observation time to the next. More specifically, a new tracking algorithm is used to compute object centric trends in the ε tot and the 11‐μmBT from two consecutive images (from the same geostationary sensor) taken at times t 1 and t 2 . In tandem with the cloud object identification technique described in the previous section, the cloud object tracking algorithm presented in this paper was specifically designed to aid in the detection of rapidly developing cloud elements in infrared satellite imagery. The object tracking methodology is not meant to be a general feature tracking solution. Attempts to utilize an existing (and available) object identification and tracking algorithm (Lakshmanan et al., 2007, 2009) were unsuccessful as the methodologies utilized in the Warning Decision Support System‐Integrated Information were unable to extract objects in a manner that consistently identified localized cloud elements within midlevel and high‐level overcast backgrounds, while being sufficiently computationally efficient for real‐time applications. x, y pixel indices). Given the variety of geostationary satellite scan strategies that are possible, data from a given geostationary satellite are mapped to a constant satellite projection in order to utilize a constant image coordinate system. Information on map projections for geostationary satellite data can be found at t 2 , the cloud object tracking algorithm consists of five major components: Define a search box and gather object properties. Compute the tracking metrics. Determine the best match‐up option. Apply quality control procedure. Compute vertical growth trends and anomalies. Object tracking logistics are simplified if done in image coordinates (pixel indices). Given the variety of geostationary satellite scan strategies that are possible, data from a given geostationary satellite are mapped to a constant satellite projection in order to utilize a constant image coordinate system. Information on map projections for geostationary satellite data can be found at http://www.cgms‐info.org/documents/cgms‐lrit‐hrit‐global‐specification‐%28v2‐8‐of‐30‐oct‐2013%29.pdf . For a given object, at time, the cloud object tracking algorithm consists of five major components: For the object of interest, at time t 2 , there are three possible outcomes of the tracking procedure: (1) The object is discarded due to lack of possible vertical growth and/or lack of proximity to a volcano, (2) matching object(s) at time t 1 is(are) identified, or (3) the object is classified as a new feature (e.g., determined not to be present in the t 1 image). Each of the five major components of the tracking algorithm is described in the subsequent sections. 3.1 Search Box and Object Properties In order to efficiently determine if a given object at time t 2 (hereafter referred to as the object of interest) is present in the t 1 image, a search box is defined in image coordinates. The size of the search box is a function of the difference between the two observation times, the horizontal resolution of the satellite data, and the attributes of the object of interest. A more detailed description of the search box can be found in Appendix A. The same search box, defined by a range of x and y pixel indices, is applied to both image times, and relevant properties of all objects, which are at least partially located in the search box, are gathered. For each of the two images, a list of objects within the search box is constructed. Each entry in the list contains the object ID, maximum value of ε tot , object area, and x, y centroids. The t 1 and t 2 lists are then utilized to construct various metrics for attempting to match the object of interest with an object or objects in the first image. However, prior to proceeding further, the algorithm first determines if the object of interest is worth tracking. For efficiency purposes, the object of interest is discarded if any of the following conditions are met: (1) all pixels are more than 200 km from the closest volcano, (2) if the maximum ε tot in the object of interest is smaller than the maximum value within the first image using the same pixels associated with the object of interest (emergent volcanic clouds should exhibit growth over an approximately fixed location), or (3) the object is likely part of a thin cirrus cloud (see Appendix B). 3.2 Object Tracking Metrics t 2 image, was deemed worthy of tracking, a cost function‐based object tracking procedure is applied. Within the search box, four tracking metrics are used to compute a cost function for each possible interimage (t 1 and t 2 ) object pair. For example, the i th t 1 object will be paired with each individual t 2 object (t 2 objects are denoted by j). This procedure is repeated for all t 1 objects in the search box to form a 2‐D matrix of cost function values. The cost function formulation for a single interimage object pair, C i,j , is shown in equation (1) If the object of interest, in theimage, was deemed worthy of tracking, a cost function‐based object tracking procedure is applied. Within the search box, four tracking metrics are used to compute a cost function for each possible interimage (and) object pair. For example, theobject will be paired with each individualobject (objects are denoted by). This procedure is repeated for allobjects in the search box to form a 2‐D matrix of cost function values. The cost function formulation for a single interimage object pair,, is shown in equation 1 C i,j are indicative of a higher‐quality interimage object match. In equation c 1 , c 2 , c 3 , and c 4 represent the four tracking metrics for a given interimage object pair i, j. The first two tracking metrics are related to object position, and the last two are related to physical properties. A combination of geometric and physical metrics is needed to ensure quality tracking for fast‐moving and/or rapidly changing objects. The first tracking metric (equation c 1 , is the fraction of pixels, in the smaller of the two objects comprising the pair (i, j) that spatially overlap (pixel counts are denoted by N). Depending on the object size and speed, spatial overlap can be strongly indicative of image‐to‐image coherence for a particular feature. However, the object tracking methodology does not require spatial overlap, nor is spatial overlap weighted any higher than the other metrics. Thus, spatial overlap does not automatically result in a valid object match. (2) Smaller values ofare indicative of a higher‐quality interimage object match. In equation 1 , andrepresent the four tracking metrics for a given interimage object pair. The first two tracking metrics are related to object position, and the last two are related to physical properties. A combination of geometric and physical metrics is needed to ensure quality tracking for fast‐moving and/or rapidly changing objects. The first tracking metric (equation 2 ), denoted by, is the fraction of pixels, in the smaller of the two objects comprising the pair () that spatially overlap (pixel counts are denoted by). Depending on the object size and speed, spatial overlap can be strongly indicative of image‐to‐image coherence for a particular feature. However, the object tracking methodology does not require spatial overlap, nor is spatial overlap weighted any higher than the other metrics. Thus, spatial overlap does not automatically result in a valid object match. c 2 ) is the normalized distance between the centroids (given by x cent and y cent ) of the object pair, where the minimum possible value of 0.0 indicates that the interimage object pair is the closest centroid‐to‐centroid distance possible given all of the possible object pairings in the search box (equation (3) (4) The second positional tracking metric () is the normalized distance between the centroids (given byand) of the object pair, where the minimum possible value of 0.0 indicates that the interimage object pair is the closest centroid‐to‐centroid distance possible given all of the possible object pairings in the search box (equation 3 ). The normalized distances are computed using the Euclidean formulation (equation 4 ). c 3 and c 4 are shown in equations ε tot and total area (A) of each object for a given interimage object pair (i, j). As equation ε tot of the t 2 object (ε tot_max,j ) and the maximum ε tot of the t 1 object (ε tot_max,i ) is normalized by whichever is larger. The relative difference in area is computed in an analogous manner (equation t 1 and t 2 . This concept will be elaborated on in the next section. (5) (6) The tracking metrics denoted asandare shown in equations 5 and 6 , respectively. Metrics 3 and 4 are aimed at quantifying the relative difference in the maximumand total area () of each object for a given interimage object pair (). As equation 5 shows, the difference between the maximumof theobject () and the maximumof theobject () is normalized by whichever is larger. The relative difference in area is computed in an analogous manner (equation 6 ). Even though the end goal is to identify rapidly developing clouds, these metrics help constrain the object tracking through influence on interimage object pairs that do have fairly similar properties atand. This concept will be elaborated on in the next section. Ideally, the cloud object tracking would consist of matching a single t 2 object with a single t 1 object. Unfortunately, that is not always the case as cloud objects frequently split and merge. For efficiency purposes, the tracking algorithm restricts splits to t 1 objects that have a ε tot_max > 0.20 and A > 200 km2 or ε tot_max > 0.30 and A > 100 km2 (the same criteria are applied to the t 2 image to similarly restrict merges). A split is defined as a single t 1 object that is associated with two or more t 2 objects. Splitting is accounted for in a given cost function entry, C i,j , when the ith t 1 object has an area that is at least twice as large as the jth t 2 object, 0.5ε tot_max,i < ε tot_max,j < 1.5ε tot_max,i (hereafter referred to as the ε tot_max split criterion), and the rectangles formed from the minimum and maximum x/y image coordinates of the interobject pair overlap (hereafter referred to as the rectangle criterion). If the splitting criteria are met, all t 2 objects that are possibly associated with the split are aggregated when constructing the C i.j cost function entry. The following criteria are used to determine which t 2 objects to aggregate: the t 2 object has an area less than the ith t 1 object, and the t 2 object satisfies the ε tot_max split and rectangle criteria. A merge is defined as multiple t 1 objects that are associated with a single t 2 object. Merging is accounted for in a given cost function entry, C i,j , when the jth t 2 object has an area that is at least twice as large as the ith t 1 object, 0.5ε tot_max,j < ε tot_max,i < 1.5ε tot_max,j (hereafter referred to as the ε tot_max merge criterion), and the rectangle criterion is satisfied. If the merging criteria are met, all t 1 objects that are possibly associated with the merge are aggregated when constructing the C i.j cost function entry. The following criteria are used to determine which t 1 objects to aggregate: the t 1 object has an area less than the jth t 2 object, and the t 1 object satisfies the ε tot_max merge and rectangle criteria. 3.3 Object Tracking Decision Logic Each value within the cost function matrix (C) corresponds to a pairing of the jth object (or jth object plus aggregate objects if split criteria are met) in the t 2 image search window and the ith object (or ith object plus aggregate objects if merge criteria are met) within the search window of the t 1 image. The cost function is analyzed with the goal of determining if the object of interest in the t 2 image can be matched with a feature in the t 1 image (with reasonable confidence) or if the object of interest is a new feature. If there are no objects within the search window of the t 1 image, the object of interest is simply classified as a new feature. If there are potential matching features in the t 1 image, then the cost function matrix is used to determine if the object of interest is a new feature or can be matched with a t 1 feature with low or high confidence. The cost function analysis begins by identifying the t 1 feature that minimizes the cost function for the object of interest (termed the primary match). Next, the t 2 feature that minimizes the cost function for the t 1 feature in the primary match is identified and designated as the secondary match. When identifying the primary and secondary matches, only interobject pairs with a normalized distance (equation 4) less than or equal to a threshold value are considered. The normalized distance threshold is determined by counting the number of possible matches between the object of interest in the t 2 image and all features in the t 1 image search window with a normalized distance less than or equal to the following values: 0.01, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, and 0.50. The threshold is taken to be the minimum value that results in at least five possible interobject pairings. If less than five interobject pairings were found after assessing all possible threshold values, a threshold of 0.50 is used. The normalized distance constraint helps prevent interobject pairs, which are likely not spatially related, from being selected in the event the cloud properties (A and ε tot_max ) are very similar. Utilizing the primary and secondary pairing information, the object of interest is classified in the following manner. If the primary match is associated with splitting or merging, the primary match is accepted and deemed to be of high quality. If the primary match is not associated with splitting or merging, the primary match is accepted, and deemed to be high quality, if the cost associated with the primary match is less than or equal to the cost function of the secondary match and the cost function of the primary match is less than a threshold. Otherwise, the primary feature pairing is classified as a low‐quality match. The cost function threshold depends on the time difference (dt) between the interimage features that compose the primary pairing and whether or not the interimage features spatially overlap. If the interimage features at least partially overlap (do not overlap), a cost function threshold of 1.25 (1.10), 1.30 (1.18), or 1.35 (1.25) is used if the dt is <18 min, 18 ≤ dt < 35, or dt ≥ 35, respectively. The cost function thresholds were determined heuristically via extensive testing of the tracking algorithm using several geostationary satellites with varying spatial and temporal capabilities. The thresholds with overlap are larger since spatial overlap adds confidence that the features are correctly paired. The 12 June 2011 eruption of Nabro (Global Volcanism Program, 2011) is used to illustrate the object identification and tracking procedure (Figure 2). In Figure 2, the ε tot_max of each object is overlaid on consecutive SEVIRI 11‐μm BT images, with starting times of 20:30 and 20:45 UTC (only a small, magnified, portion of the full disk image is shown to allow for detailed analysis centered on Nabro). The object ε tot_max , at each time, is assigned to all pixels that comprise a given object. The boundaries of the objects are evident by the presence of the underlying gray scale 11‐μm BT image relative to the colorized ε tot_max values. Objects that were deemed worthy of tracking (see section 3.1) are labeled with a unique identification number. For clarity, only tracking‐worthy objects that are composed of three or more pixels at 20:45 UTC are labeled. Using the automated tracking procedure, a feature (more than one object if a split was detected) with a given identification number in the t 2 image (Figure 2b) was matched with a feature (more than one object if a merge was detected) with the same identification number in the t 1 image (Figure 2a). Inspection of Figure 2 reveals that the feature‐tracking algorithm generally produces feature pairings that are consistent with human assessment. Further, Figure 3 (a–d) shows several object tracking parameters overlaid on the t 2 SEVIRI image (20:45 UTC). The object tracking classification flag (Figure 3a) shows that many of the objects are not tracked due to being classified as thin cirrus, not having any possibility of vertical development, or being more 200 km from the nearest volcano. The labeled objects in the 20:45 UTC image are generally matched with features in the 20:30 UTC image with high confidence, including two correctly identified merges (Features 3 and 5). Given the relatively small time interval (15 min) between the images, most of the feature pairings spatially overlap (Figure 3b) and have centroids that are closer than alternate object pairings (Figure 3c). The two low confidence pairings (Features 6 and 12) are associated with small transient objects with large tracking cost function values (Figure 3d). Figure 2 Open in figure viewer PowerPoint The maximum top of troposphere emissivity (ε tot_max ) associated with each objectively defined cloud object is overlaid on gray scale SEVIRI 11‐μm brightness temperature imagery from (top) 20:30 and (bottom) 20:45 UTC on 12 June 2011. Interimage cloud object pairings, derived from the automated tracking algorithm, are labeled with a unique identification number if the cloud object is composed of 3 or more pixels at 20:45 UTC. Figure 3 Open in figure viewer PowerPoint Several cloud object tracking parameters are overlaid on gray scale SEVIRI 11‐μm brightness temperature (BT) imagery from 20:45 UTC on 12 June 2011. The tracking classification (low confidence 1‐to‐1 match, high confidence 1‐to‐1 match, high confidence split, high confidence merge, new object, untracked cirrus, untracked objects due to lack of possible growth, untracked objects due to large distance from nearest volcano, and objects discarded due to failed quality control) is shown in panel a. The interimage object overlap fraction, distance ratio, and tracking cost function are shown in panels b, c, and d, respectively. The change in the minimum cloud object 11‐μm BT and the vertical growth anomaly are shown in panels e and f. The identification number of each tracked feature is also shown on each panel. 3.4 Cloud Trends The volcanic eruption detection algorithm is largely predicated on identifying clouds that exhibit anomalous vertical growth from one image to the next. Vertical growth, or lack thereof, is assessed using the difference between the ε tot_max of the t 2 feature and the corresponding representation of the t 2 feature in the t 1 image. The difference between the minimum BT (BT min ) of the t 2 feature and the corresponding representation in the t 1 image is also computed. Given that the cloud object tracking methodology is restricted to dt < 60 min, changes in ε tot_max (dε tot ) and BT min (dBT) will largely be associated with changes in cloud height and/or opacity, as opposed to changes in temperature within a cloud layer with approximately constant vertical boundaries. If the t 2 feature of interest was confidently traced back to a t 1 feature, using the methodology described in section 3.3, dε tot and dBT are simply computed using the difference between the matching t 2 and t 1 features, such that positive (negative) values of dε tot (dBT) generally correspond to vertical growth. If the t 2 feature was classified as a new feature, the t 1 ε tot_max (BT min ) is given by the maximum (minimum) value within the t 1 pixels that correspond to the location of the t 2 object (termed the object footprint). In addition, the object footprint approach is also applied if the t 1 and t 2 features were matched with low confidence and the difference between the ε tot_max in the t 2 feature of interest and the maximum ε tot within the search window of the t 1 image is greater than 0.05, which is consistent with a t 2 feature that has undergone vertical development. Given the goal of identifying clouds associated with explosive volcanic eruptions, the object footprint approach to computing dε tot and dBT is applied to low confidence matches associated with cloud vertical growth, to ensure an accurate result if the cloud formation is directly linked to volcanic activity (e.g., cloud formation centered on a volcanic vent). If the t 2 feature of interest has a low confidence t 1 pairing and does not have a ε tot_max that is at least 0.05 greater than the ε tot_max in the t 1 search window, the t 1 ε tot_max and BT min is taken to be the maximum and minimum value, respectively, within the t 1 search window. Figure 3e shows the dBT for the Nabro volcanic cloud and surrounding meteorological clouds. The Nabro cloud has a dBT of −45.3 K, which is emblematic of very strong vertical growth in the 15‐min time window. None of the surrounding meteorological clouds had a comparable decrease in BT min . 3.5 Cloud Trend Quality Control Cirrus clouds can be spatially expansive, fast moving, and/or transient, making them difficult to track. Thus, a simple quality control procedure is applied with the aim of identifying bad or broken cloud tracks, which are most common when cirrus cloud elements are being tracked. Unfortunately, not all cirrus can be effectively screened out using the method described in Appendix B, but the screening procedure is useful. For each cloud object that exhibits vertical growth between two successive satellite images, the mean ε tot of the pixels that reside just outside of the object of interest is computed for both images. Outside pixels are defined as those that are not part of the object of interest but are within a rectangular region defined by the minimum and maximum x and y coordinates of the object of interest at a given time, with an additional 2 pixel buffer applied. The rectangular region within the first image will differ from the second image if the spatial footprint of the object differs between the two images, which it usually does. The mean ε tot of the outside pixels from the first image is compared to the second image, with the expectation that the ε tot of the pixels immediately surrounding the object of interest will only increase significantly from the first image to the second if the object track is in error in the presence of a complex and partially transient object field with large spatial gradients in ε tot . Regions with fast‐moving and rapidly evolving cirrus, such as near jet streams, are an example of a highly complex object field. The object tracking fails the quality control if the mean ε tot of the outside pixels increases by 95%, 45%, 25%, or 10%, depending on if the mean outside pixel ε tot of the second image is greater than 0.07, 0.10, 0.30, or 0.50, respectively. The quality control thresholds were determined heuristically upon extensive analysis of cirrus‐related tracking errors identified through near‐real‐time processing of several satellite data streams. This simple quality control procedure was found to eliminate about half of the cirrus‐related tracking errors, with limited risk of filtering out genuine cloud vertical growth.

4 Meteorological Cloud Analysis rings in Figure GOES‐13 Imager (depicted as GOES‐East in Figure 4)—One complete day from every month of 2012 (image interval varies from 5 to 30 min) was utilized. Several days were selected to coincide with intense meteorological convection as indicated by the National Oceanic and Atmospheric Administration (NOAA) Storm Prediction Center storm report database (http://www.spc.noaa.gov/). MSG‐2 SEVIRI full disk—One complete day from every month of 2012 (image interval is 15 min) was utilized. Several days were selected to coincide with intense meteorological convection as indicated by the European Severe Storms Laboratory storm report database (https://www.essl.org). MSG‐1 SEVIRI European mesoscale sector—Four complete days, with convectively active weather, in June 2009 (image interval is 5 min) were chosen. MTSAT‐2 Imager—One complete day from each month between January and October 2012 (image interval varies from 15 to 60 min) was utilized. November and December were excluded due to data availability issues. Himawari‐8 AHI full disk—One complete day from each month between June 2015 and May 2016 was chosen for the AHI full disk (image interval is 10 min). Himawari‐8 AHI Japan region sector—One complete day from each week between 14 April 2016 and 10 June 2016 was chosen for the AHI Japan sector (image interval is 2.5 min). In order to detect a wide range of explosion volcanic eruptions, the cloud vertical growth rates must be compared to cloud vertical growth associated with meteorological processes. However, satellite‐inferred cloud vertical growth depends not only on the vigor of the upward motion but also on observational factors. The primary observational factors are the time interval between the successive images, the area of the satellite pixels, and the initial cloud state. Taking into account the primary observational factors, a database of cloud vertical growth associated with meteorological processes was generated. The database is composed of a large volume of data from multiple satellites (see Table 1 for acronym definitions) with significantly varying spatial resolution and image refresh intervals. The database was constructed to be as geographically and seasonally complete as possible. The geographic extent is depicted by the geostationary coveragein Figure 4 (the volcano names in Figure 4 are unrelated to the meteorological database component of this study). The database entries were chosen such that no volcanic ash advisories for clouds above 20,000 feet (~6 km) were in effect and the Global Volcanism Program ( http://volcano.si.edu/ ) also indicated no high‐level volcanic clouds. The following data sets were chosen for inclusion into the database: Table 1. The Geostationary Satellite Imaging Radiometers Utilized in This Study Sensor acronym Acronym meaning GOES Imager Geostationary Operational Environmental Satellite (GOES) Imager GOES‐R ABI Next Generation Geostationary Operational Environmental Satellite (GOES‐R) Advanced Baseline Imager (ABI) Himawari‐8/Himawari‐9 AHI Himawari‐8/9 Advanced Himawari Imager (AHI) MTSAT Imager Multifunctional Transport SATellites (MTSAT) Imager MSG SEVIRI Meteosat Second Generation (MSG) Spinning Enhanced Visible and Infrared Imager Figure 4 Open in figure viewer PowerPoint The colorized rings illustrate the approximate coverage area of the geostationary satellites utilized in this study.R elevant volcanoes are also labeled. Satellite acronyms are defined in Table 1 Histograms of cloud vertical growth were constructed from each meteorological cloud object that exhibited vertical growth between two successive images. A total of 6,453,781 t 2 objects exhibited vertical growth between two successive images. Many histograms were constructed to account for the time interval among successive images, the pixel area, and the initial cloud state (characterized by ε tot_max at t 1 ). The binning scheme used to account for the observational factors is given in Table 2. For each growth rate distribution, the mean and standard deviation were computed. Figure 5 helps illustrate the sensitivity to the time interval between successive satellite images (dt). For a given pixel area and initial cloud state, the smaller the time interval between images, the greater constraint on the actual time interval over which the vertical development occurred. The time sampling effect results in a narrowing of the dBT distribution as the sampling interval decreases. Thus, a given strong decrease in the BT min of a cloud object will be increasingly anomalous with decreasing dt. This suggests that next generation satellites such as Himawari‐8 and GOES‐R, with more frequent sampling, are well suited for flagging strong vertical growth anomalies. Most vertically growing clouds are spatially inhomogeneous, with the BT min located within a localized core. Thus, the observed BT min (and ε tot_max ) is a strong function of the spatial resolution of the satellite measurements, where the satellite‐inferred vertical extent of the cloud core approaches the true vertical extent as the satellite pixel area decreases. As the pixel area increases, cumuliform clouds become increasingly spatially uniform from the satellite perspective, which results in a narrower distribution of dBT, as shown in Figure 6 (the image time interval and the initial cloud state are constrained in Figure 6). For a given viewing angle, sensors with improved spatial resolution, such as Himawari‐8 and GOES‐R, require larger decreases in BT min to produce a given growth anomaly compared to heritage sensors with poorer spatial resolution. The impacts of the initial cloud state (Figure 7) are a bit more complicated. For a given image time interval and pixel area, the width of the dBT distribution generally decreases as the initial cloud state becomes increasingly indicative of midlevel and high‐level clouds. The exception occurs when the meteorological cloud evolves from a clear‐sky or low‐level cloudy state, characterized by a small ε tot_max at t 1 .The narrower dBT distribution at small values of ε tot_max (at t 1 ) is likely a reflection of the atmosphere not always being favorable for significant convection. Table 2. The Binning Scheme Used to Account for the Impact of Observational Factors on Cloud Vertical Growth Rate Calculations Observational factor Binning scheme Comments Image interval (min) 7 bins: The binning scheme is designed to account for many different scan capabilities and strategies. 1–4, 4–7, 7–11, 11–18, 18–27, 27–35, >35 Pixel area (km2) 6 bins: The binning scheme accounts for varying pixel area associated with sensor capabilities and viewing angle. 4–6, 6–10, 10–20, 20–30, 30–50, >50 Background conditions (given by the maximum ε tot of the matching object in the first image of the image pair) 19 bins: This binning scheme accounts for the initial cloud state (or lack there of) of a cloud object that exhibits vertical growth between two successive satellite images. 0.0–0.90 in increments of 0.05, and a final bin of >0.90 Figure 5 Open in figure viewer PowerPoint Distribution of meteorological cloud vertical growth rates (dBT(11 μm)), defined as the decrease in the minimum cloud object 11‐μm brightness temperature (BT) between two successive images, is shown for measurements taken at 2.5‐ (panel a) and 10‐min (panel b) intervals. The dashed blue lines denote the 2, 4, 6, 8, and 10 standard deviation anomaly values. These distributions were constrained to consist of cloud objects with similar pixel area and initial cloud state. Figure 6 Open in figure viewer PowerPoint 2 are shown in panels a, b, c, and d, respectively. Same as Figure 5 , except that the cloud vertical growth rate distributions are shown as a function of pixel area while constraining the time interval between successive images and the initial cloud state. Cloud objects with a mean pixel area of 4–6, 6–10, 10–20, and 20–30 kmare shown in panels a, b, c, and d, respectively. Figure 7 Open in figure viewer PowerPoint ε tot , at the first image time in an image pair, of 0.00–0.05, 0.50–0.55, 0.70–0.75, and >0.90 are shown in panels a, b, c, and d, respectively. Same as Figure 5 , except that the cloud vertical growth rate distributions are shown as a function of the initial cloud state while constraining the time interval between successive images and the pixel area. Cloud objects with a maximum, at the first image time in an image pair, of 0.00–0.05, 0.50–0.55, 0.70–0.75, and >0.90 are shown in panels a, b, c, and d, respectively. The meteorological cloud database, which is function of the dt, pixel area, and initial cloud state, is used to compute the standardized vertical growth rate anomaly (subtract the mean and divide by the standard deviation) of cloud objects that exhibit vertical growth between two successive images. Anomalies are computed for both dBT (the dBT anomaly is multiplied by −1) and dε tot , and the larger of the two anomalies is taken. The dBT growth metric ensures that undercooled (colder than the tropopause) clouds, often associated with stratospheric tops, are properly characterized, while the dε tot metric is better at characterizing growth in atmospheres with a relatively small tropospheric vertical temperature gradient. Example standardized vertical growth rate anomalies, hereafter referred to as z‐scores, are shown in Figure 3f. The Nabro eruptive cloud has a z‐score of about 8, which is significantly greater than any of the surrounding meteorological clouds. Many more volcanic eruptions were analyzed, as described in the next section.

5 Volcanic Eruption Analysis The CGA method, which automatically identifies, tracks, and characterizes (quantifies cloud vertical growth relative to a database of meteorological clouds) cloud objects, was applied to more than 70 different explosive volcanic events from 30 volcanoes with comprehensive global sampling (see volcano labels in Figure 4 and Table 3 through Table 6). Results from 12 different geostationary satellites, with varying spatial coverage and capabilities, ranging from heritage to next generation, were included in the analysis. Prior to analyzing the results, the performance of the cloud object identification and tracking algorithm was assessed by comparing the autogenerated results to a manual analysis of images similar to those shown in Figures 1–3. The automated algorithm produced cloud vertical growth metrics (dBT and dε tot ) that were predominantly consistent with the manual assessment. Disagreements with the manual analysis were found in only 10 cases, as indicated by table entries with two sets of growth rate anomaly statistics (see Table 3 through Table 6). In 9 of those 10 cases, the volcanic object of interest was matched with an unrelated meteorological cloud object. The mismatch errors result in the cloud growth anomaly being underestimated or overestimated, usually by no more than 40%. As described in the section 6, the practical impact of the mismatch errors is generally small. Most (eight of nine) of the mismatch errors were associated with legacy satellites with less frequent imaging, such as GOES‐8/GOES‐9/GOES‐10/GOES‐11/GOES‐12/GOES‐13 and MTSAT‐1R/2. For instance, the growth rate anomaly of the 31 July 2015 eruption of Manam (Papua New Guinea, PNG), as observed in 60‐min intervals by MTSAT‐2, was underestimated by 41% due to the volcanic cloud being matched with an unrelated convective meteorological cloud. Given the 60‐min image interval, such an object tracking error, in a convectively active region, is not surprising. The same eruption was accurately characterized with 10‐min Himawari‐8 imagery (see Table 3 through Table 6), which helps illustrate the value of more frequent imaging. The 4 January 2016 (13:00 UTC) Soputan (Indonesia) case, captured by Himawari‐8, was the only mismatch error associated with next generation satellite capabilities (see Figure 8). The automated tracking algorithm matched the Soputan cloud (Feature 8 on the right side panels of Figure 8) with a meteorological cloud west of Soputan in the 12:50 UTC Himawari‐8 image. The resulting z‐score was 15.9. In contrast, the z‐score was 35.7 if the Soputan cloud is assumed to have grown in place. Animations of Himawari‐8 imagery show that the meteorological cloud matched with the Soputan cloud by the automated algorithm may have advected into the pixels occupied by the Soputan cloud, although this assessment is inconclusive. The Soputan case, though, does illustrate the potential for false amplification of vertical growth due to cirrus advection over developing cumulus clouds. The false amplification side effect is, at least partially, mitigated if the tracking algorithm matches the vertically developing cloud with the correct cirrus cloud in the first image of the image pair. Nevertheless, the z‐score of 15.9, retrieved by the automated algorithm, is sufficiently large to warrant attention, even if underestimated. Thus, from a practical standpoint, the tracking algorithm correctly highlights the Soputan cloud as a feature of interest, even if the cloud tracking process was imperfect. Table 3. A Summary of Cloud Vertical Growth Associated With Explosive Volcanic Events Captured by GOES Satellites Volcano Country Date mm/dd/yy UTC Time UTC dt min dBT °C dε tot z GOES‐8 Reventador0 Ecuador 11/03/02 12:45 30.0 −28.7 0.31 4.4 GOES‐9 Manam0 PNG 10/24/04 01:25 60.0 −59.5 0.53 7.8 Manam0 PNG 01/27/05 14:02 36.3 −59.1 0.53 6.9 −43.4 0.36 5.3 Anatahan0 USA 04/05/05 16:13 48.0 −37.9 0.44 6.1 GOES‐10 Chaitén0 Chile 05/02/08 12:28 12.9 −3.9 0.04 −0.1 Chaitén0 Chile 05/02/08 12:45 17.0 −14.8 0.24 4.7 Chaitén1 Chile 05/02/08 15:45 17.1 −39.9 0.49 10.4 Chaitén1 Chile 05/06/08 12:28 13.0 −35.8 0.56 12.9 −41.3 0.68 12.0 Chaitén2 Chile 05/06/08 18:45 17.1 −28.1 0.39 8.9 −22.2 0.32 6.4 Chaitén0 Chile 05/07/08 00:15 17.5 −12.5 0.18 3.5 Chaitén0 Chile 05/07/08 02:45 29.5 −23.2 0.45 4.8 Chaitén0 Chile 05/07/08 20:45 29.5 −20.8 0.32 3.4 Chaitén0 Chile 05/08/08 03:45 17.0 −17.9 0.30 4.8 Galeras0 Colombia 03/13/09 21:15 19.7 −15.8 0.17 1.7 GOES‐11 Okmok1 USA 07/12/08 20:00 15.4 −35.0 0.49 17.6 Kasatochi2 USA 08/07/08 22:30 14.2 −22.7 0.29 7.4 Redoubt1 USA 04/04/09 14:15 14.6 −26.8 0.45 16.0 GOES‐12 Santa Ana1 El Salvador 10/01/05 14:45 32.0 −67.4 0.70 10.6 Galeras0 Colombia 02/20/09 12:15 27.6 −32.7 0.41 7.0 Puyehue‐Cordón2 Caulle Chile 06/04/11 18:28 12.9 −30.4 0.46 8.2 −24.5 0.38 6.0 GOES‐13 Soufrière Hills1 Montserrat 02/11/10 17:15 30.0 −64.4 0.61 18.1 Popocatépetl2 Mexico 04/14/12 03:32 14.3 −32.4 0.38 7.9 Nevado del Ruiz1 Colombia 05/29/12 08:15 30.0 −42.5 0.43 5.8 −44.5 0.51 6.2 Tungurahua1 Ecuador 07/14/13 11:45 32.5 −64.6 0.81 23.5 Calbuco1 Chile 04/22/15 21:38 33.8 −60.4 0.84 17.3 Calbuco1 Chile 04/23/15 04:08 30.0 −59.7 0.84 14.5 Wolf1 Ecuador 05/25/15 07:15 30.0 −52.2 0.63 18.8 Cotopaxi1 Ecuador 08/14/15 09:15 27.3 −44.8 0.51 6.2 −49.6 0.62 9.7 Table 4. Same as Table , Except for Volcanic Events Captured by the MTSAT Satellites Volcano Country Date mm/dd/yy UTC Time UTC dt min dBT °C dε tot z MTSAT‐1R Sarychev Peak1 Russia 06/12/09 01:57 27.0 −39.2 0.56 10.7 Sarychev Peak1 Russia 06/13/09 21:30 33.0 −47.0 −50.2 0.70 0.76 18.7 14.2 Sarychev Peak1 Russia 06/14/09 18:57 27.0 −44.1 0.69 18.3 Sarigan0 USA 05/29/10 12:30 60.0 −56.7 0.51 7.5 Kelut1 Indonesia 02/13/14 16:19 10.0 −53.2 0.28 40.2 Tongariro0 NZ 08/06/12 12:01 29.0 −26.0 0.38 9.8 MTSAT‐2 Kelut0 Indonesia 02/13/14 16:32 60.0 −37.0 0.31 6.9 Sangeang Api1 Indonesia 05/30/14 08:32 60.0 −96.2 0.96 25.8 Sheveluch0 Russia 03/25/15 22:32 31.0 −21.0 0.27 3.6 Manam0 PNG 07/31/15 01:32 60.0 −35.6 −65.9 0.27 0.65 6.6 11.2 Table 5. Same as Table , Except for Volcanic Events Captured by the MSG Satellites Volcano Country Date mm/dd/yy UTC Time UTC dt min dBT °C dε tot z Meteosat‐8 (MSG‐1) Karthala2 Comoros 11/24/05 18:15 15.0 −25.6 0.25 6.7 Grímsvötn 1 Iceland 05/21/11 19:19 5.0 −24.6 0.35 21.3 Meteosat‐9 (MSG‐2) Soufrière Hills1 Montserrat 02/11/10 17:15 15.0 −52.2 0.45 15.1 Grímsvötn1 Iceland 05/21/11 19:15 15.0 −26.4 0.36 12.7 Nabro1 Eritrea 06/12/11 20:45 15.0 −45.3 0.45 8.1 Meteosat‐10 (MSG‐3) Etna1 Italy 11/23/13 09:45 15.0 −24.0 0.42 11.3 Etna2 Italy 11/23/13 10:00 15.0 −30.6 0.43 7.4 Etna2 Italy 12/03/15 02:15 15.0 −17.8 0.29 7.4 Etna1 Italy 12/03/15 02:30 15.0 −41.9 0.60 10.9 Etna0 Italy 12/04/15 09:00 15.0 NA NA NA −18.9 0.28 6.2 Etna2 Italy 12/04/15 09:15 15.0 −35.5 0.53 8.7 Etna2 Italy 12/04/15 20:45 15.0 −25.4 0.41 6.9 Etna2 Italy 12/05/15 15:00 15.0 −26.8 0.38 6.0 Table 6. Same as Table , Except for Volcanic Events Captured by the Himawari‐8 Satellite Volcano Country Date mm/dd/yy UTC Time UTC dt min dBT °C dε tot z Himawari‐8 Sheveluch1 Russia 03/25/15 22:10 10.0 −32.0 0.49 10.4 Kuchinoerabujima2 Japan 05/29/15 01:10 10.0 −16.3 0.22 8.6 Manam1 PNG 07/31/15 01:30 10.0 −62.1 0.51 10.5 Soputan1 Indonesia 01/04/16 13:00 10.0 −51.5 −63.7 0.56 0.78 15.9 35.7 Soputan1 Indonesia 01/04/16 22:40 10.0 −20.0 0.30 10.9 Soputan1 Indonesia 01/04/16 22:50 10.0 −36.0 0.42 8.5 Soputan1 Indonesia 01/05/16 06:40 10.0 −39.4 0.43 9.3 Soputan1 Indonesia 01/05/16 15:20 10.0 −23.8 0.35 13.2 Bogoslof2 USA 12/21/16 00:50 10.0 −13.1 0.09 5.2 Bogoslof1 USA 12/22/16 01:30 10.0 −19.3 0.25 11.2 Bogoslof1 USA 12/26/16 23:30 10.0 −31.0 0.36 22.1 Bogoslof0 USA 12/30/16 08:50 10.0 −21.0 0.14 3.9 Bogoslof1 USA 01/04/17 06:30 10.0 −38.4 0.40 25.0 Bogoslof2 USA 01/05/17 22:30 10.0 −13.1 0.12 6.0 Bogoslof2 USA 01/09/17 08:10 10.0 −14.6 0.20 7.7 Bogoslof2 USA 01/18/17 22:30 10.0 −15.3 0.21 9.1 Bogoslof0 USA 01/20/17 22:20 10.0 −6.7 0.08 3.0 Bogoslof0 USA 01/22/17 23:10 10.0 −7.9 0.12 4.3 Bogoslof1 USA 01/24/17 14:00 10.0 −20.9 0.20 11.7 Bogoslof1 USA 01/26/17 16:00 10.0 −20.9 0.24 10.8 Bogoslof2 USA 01/27/17 17:40 10.0 −16.1 0.16 7.0 Bogoslof0 USA 01/31/17 06:50 10.0 −1.9 0.03 0.7 Bogoslof0 USA 02/18/17 00:40 10.0 −10.0 0.08 4.4 Bogoslof1 USA 02/18/17 14:10 10.0 −18.0 0.27 12.2 Bogoslof1 USA 02/20/17 02:20 10.0 −25.5 0.24 14.3 Bogoslof0 USA 03/08/17 08:20 10.0 −12.3 0.07 4.9 Bogoslof0 USA 03/13/17 11:40 10.0 −9.0 0.07 4.7 Sheveluch2 Russia 09/07/17 12:50 2.5 −11.0 0.12 5.2 Sheveluch1 Russia 09/07/17 13:00 10.0 −25.0 0.26 11.3 Figure 8 Open in figure viewer PowerPoint (first row) Himawari‐8 false color imagery, (second row) 10.4‐μm brightness temperature (BT), and (third row) 10.4‐μm top of troposphere emissivity (ε tot ) are shown for two consecutive image times. The first column is from the 12:50 UTC Himawari‐8 AHI scan on 4 January 2016. The second column is from the 13:00 UTC Himawari‐8 AHI scan on 4 January 2016. The perimeter of automatically defined cloud objects is overlaid on each image in black, and object tracking results can be inferred by matching the integer labels from the 12:50 UTC image and the 13:00 UTC image. The other major algorithm issue was associated with the 4 December 2015 eruption of Etna (Italy), which produced a cloud that first became visible in MSG‐3 imagery at 9:00 UTC. Manual analysis of the imagery indicated that the 9:00 UTC Etna cloud had a z‐score of 6.2. The 9:00 UTC Etna cloud was classified as cirrus in the preprocessing step of the object tracking algorithm and, hence, was not tracked, and the vertical growth was not characterized. The cirrus classification, while not practically desirable, is accurate since multispectral imagery reveals the presence of widespread cirrus clouds above Etna and the surrounding areas. The Etna cloud initially appeared to develop underneath the cirrus shield. The Etna cloud vertical growth, however, was effectively computed in the 9:15 UTC MSG‐3 image (vertical growth that occurred between the 9:00 and 9:15 UTC images was properly captured). Thus, sometimes, event detection with the GCA method will be slightly delayed when volcanic clouds, from explosive eruptions, initially develop underneath meteorological clouds. As will be discussed in section 6, the cloud object tracking algorithm also does occasionally result in false vertical growth being detected. The CGA method is designed to aid in timely eruption detection; it is not intended to be a tool for classifying eruptions. No effort is made to relate cloud growth anomalies to metrics such as the Volcanic Explosivity Index (Newhall & Self, 1982), as the cloud growth anomalies are sensitive to nonvolcanic factors such as the background atmospheric state and the time of the satellite observations relative to the evolution of the volcanic cloud. The following analysis is designed to illustrate how clouds produced by explosive volcanic eruptions evolve relative to meteorological clouds and assess the impact of satellite sensor capabilities on quantifying cloud vertical growth. Figure 9a shows the distribution of z‐scores for 79 volcanic clouds processed using the CGA method. Some eruptions were observed with multiple sensors and/or produced clouds that grew significantly over more than one image pair (see Table 3 through Table 6). Out of the 79 total cases, the z‐score was found to be 5 or greater in 72 cases and 10 or greater in 22 cases. Thus, volcanic clouds produced by explosive eruptive events (maximum cloud altitude in the upper troposphere or higher) systematically develop in an anomalous manner compared to meteorological clouds. All of the cases with z‐scores of 15 or greater are shown in Figure 9b. The top 15 eruption‐related anomalies, shown in Figure 9b, are about evenly split between tropical and nontropical volcanoes. The sampling is insufficient for determining if tropical eruptions are more or less likely to produce large anomalies compared to nontropical eruptions. Figure 9 Open in figure viewer PowerPoint The distribution of cloud vertical growth rate anomalies, expressed as the number of standard deviations above the mean growth rate derived from a meteorological cloud data set, is shown in panel a for 79 volcanic cases. The cumulative frequency (black line) is also shown relative to the y axis on the right. A ranking of the top 15 volcanic vertical growth rate anomalies (red bars) is shown in panel b, along side growth rate anomalies associated with a supercell thunderstorm and the development of wildfire‐driven convective storm (blue bars). The analysis reveals that satellite capabilities can have a profound impact on the z‐scores, with the CGA method being more effective when applied to next generation satellites, such as Himawari‐8. For instance, the 25 March 2015 eruption of Sheveluch (Russia) produced a z‐score of 10.4 with Himawari‐8 10‐min imagery, while the same event had a z‐score of 3.6 when the growth is derived from MTSAT‐2 images taken 31 min apart (see Table 3 through Table 6). The larger Himawari‐8 z‐score can be mostly attributed to a combination of more frequent imaging (thereby constraining the time period of cloud growth) and the factor of 4 improvement in spatial resolution (more detailed imaging of cloud core). The Himawari‐8‐derived z‐score for the 31 July 2015 eruption of Manam (PNG) was 10.5, whereas the corresponding MTSAT‐2 z‐score was 6.6. The interval between MTSAT‐2 images was 60 min, compared to 10 min for Himawari‐8. The time interval between successive satellite images seems to be the most critical measurement attribute that impacts the CGA method. The 2014 Kelut (Indonesia) and 2011 Grímsvötn (Iceland) events illustrate the value of increased image refresh. The onset of the Kelut eruption was observed by MTSAT‐1R at 10‐min intervals and by MTSAT‐2 at 60‐min intervals. The spatial resolution of MTSAT‐1R and MTSAT‐2 at Kelut is very similar, and the spectral bands on each sensor are nearly identical. Thus, time sampling dominates the difference in measurement capabilities. When imaged every 60 min with MTSAT‐2, the z‐score was 6.9. In contrast, the 10‐min MTSAT‐1R images resulted in a z‐score of 40.2, which is the largest z‐score found to date using the CGA method. The Grímsvötn event was captured by MSG‐1 at 5‐min intervals and MSG‐2 at 15‐min intervals. Other than the difference in the time sampling, the MSG‐1 and MSG‐2 views of Grímsvötn are very similar. The z‐score associated with the 15‐min imagery was 12.7, compared to 21.3 for the 5‐min imagery. The Grímsvötn example illustrates that increased time sampling can have a significant impact, even when images are already refreshed relatively frequently. It is worth noting that the 2010 Soufrière Hills event produced similar z‐scores when viewed by GOES‐13 at 30‐min intervals and MSG‐2 at 15‐min intervals. The similarity is largely a result of the eruption occurring in a nearly cloud‐free environment combined with the volcanic cloud development being confined to the second half of the GOES‐13 30‐min sampling interval. Thus, the impact of time sampling will vary depending on the background conditions and eruption attributes. However, there is no doubt that more frequent imaging significantly improves the effectiveness of the CGA method through more accurate cloud tracking and cloud vertical development characterization. For added perspective, the vertical growth anomalies associated with volcanic eruptions were compared to known extreme cases of cloud growth unaffiliated with volcanism. Wildfire induced cumulonimbus (pyroCb) clouds are known to develop rapidly (Fromm et al., 2006, 2010). A notable pyroCb was observed by MSG‐2 on 4 August 2010 (http://oiswww.eumetsat.org/WEBOPS/iotm/iotm/20100804_pyrocb/20100804_pyrocb.html). The pyroCb was the result of a large wildfire 300 km southeast of Moscow, during a pronounced Russian heat wave. When the CGA method was applied to the pyroCb, the maximum z‐score, depicted by the blue bar labeled 17 in Figure 9b, was 4.1. About 80% of the volcanic eruptions analyzed had a larger z‐score. The 2010 eruption of Soufrière Hills was observed by MSG‐2 at roughly the same viewing angle and background cloud cover as the Russian pyroCb. In contrast, the Soufriere Hills cloud had a maximum z‐score of 15.1. A developing supercell thunderstorm, which subsequently produced an EF‐5 tornado near El Reno, OK, USA, on 24 May 2011 (Bluestein et al., 2015), was also analyzed. The GOES‐13 satellite captured the early development of the El Reno storm at 7–10‐min imaging intervals. The maximum observed z‐score, which is shown in Figure 9b, was 11. While only about 27% of the volcanic eruptions analyzed have a larger z‐score (Figure 9a), the El Reno storm was an extreme outlier compared to the typical supercell (Cintineo et al., 2013). Supercells of such intensity are fortunately rare, and the favored locations for supercell development (large plains) are generally in volcanically inactive regions. Thus, from a practical perspective, supercell development is not a common source of ambiguity for volcanic eruption identification using the CGA method. Automated eruption detection and false alarm sources are discussed in detail in the next section.

6 Alerting Tool In order to be useful for automated alerting, the CGA method must not only be capable of reliably detecting explosive volcanic events, with clouds that reach the upper troposphere or higher, but must do so with a manageable number of false alarms. The CGA method is not designed to detect all volcanic eruptions, so ideally, it should be combined with complementary detection methods (e.g., Brenot et al., 2014; Pavolonis et al., 2015a, 2015b), which is a topic of future research. In order to generate volcanic cloud alerts from the CGA products, thresholds must be used to select vertically growing cloud objects of interest that are in close proximity to known volcanic vents. The selection thresholds are given in Table 7. The selection thresholds were determined heuristically using the previously analyzed 79 volcanic events and results from continuous near‐real‐time processing that commenced in 2013. In order to account for increases in cloud tracking uncertainty as the time between successive satellite images increases, the selection criteria are a function of the time interval between images. The selection criteria also place greater weight on cloud objects that are closer to volcanic vents after a parallax correction has been applied. A combination of dBT and the z‐score is used to ensure that the computed anomalies are truly associated with noteworthy vertical growth. Finally, less stringent selection thresholds can optionally be applied to cloud objects that are near volcanoes with known ongoing activity. For real‐time application of the CGA, a list of currently active volcanoes, as given by the Global Volcanism Program, is maintained. Table 7. Criteria Used to Identify Potential Volcanic Eruptions From Cloud Vertical Growth dt (min) dBT (°C) z r (km) Rε tot 65 −80 25 75 NA 65 −70 10 75 0.01 65 −50 10 25 0.90 35 −45 10 25 0.01 35 −35 10 25 0.20 35 −40 5 25 0.75 18 −45 2 25 0.01 18 −35 10 25 0.10 18 −25 10 25 0.25 18 −20 10 25 0.50 18 −15 5 25 0.01* 18 −15 5 10 NA* 11 −30 8 25 0.01 11 −15 10 25 0.01 11 −10 5 25 0.01* 11 −10 5 10 NA* 6 −20 8 25 0.01 6 −10 10 25 0.01 The CGA algorithm and CGA‐based eruption alerting capability are part of the VOLcanic Cloud Analysis Toolkit (VOLCAT). VOLCAT is a collection of software developed by the NOAA, in partnership with the University of Wisconsin‐Madison. In addition to the CGA algorithm, VOLCAT includes the volcanic cloud‐related algorithms described in Pavolonis et al. (2013, 2015a, 2015b). Example CGA‐based alerts from VOLCAT are shown in Figure 10. The alerts are in the form of a hyperlink, distributable via e‐mail or Short Message Service. The hyperlink points to a web‐based alert report that includes information on the cloud growth anomaly, a list of most likely source volcanoes, and relevant satellite imagery. Satellite imagery, including links to animations, is included to aid in alert interpretation. The example alerts shown in Figure 10, which were generated in an unsupervised near‐real‐time manner, capture the very early evolution of clouds produced by events at Sheveluch (Russia) and Ambae (Vanuatu). Both alerts were based on 10‐min Himawari‐8 data. The VOLCAT alerts are currently distributed to expert users at Volcanic Ash Advisory Centers and Volcano Observatories in an experimental manner. Figure 10 Open in figure viewer PowerPoint Example volcanic cloud alerts that were automatically generated in near real time on (top) 11 November 2017 and (bottom) 29 October 2017 are shown. The alert in the top panel captured an eruption of Sheveluch (Russia) and the bottom panel captured an eruption of Ambae (Vanuatu). Both eruptions were detected in a timely manner using Himawari‐8 satellite data and the cloud growth anomaly detection algorithm. To gain insight into the typical performance of the alerting system, all CGA‐based VOLCAT alerts were manually analyzed for a 32‐day period from 6 October to 6 November 2017. During that period, eight confirmed CGA alerts were generated for events at Popocatépetl (Mexico), Sheveluch (Russia), and Ambae (Vanuatu). The two Popocatépetl events (7 October at 16:27 UTC and 4 November at 20:47 UTC) were detected using preliminary nonoperational GOES‐16 data with 5‐min image refresh. GOES‐16, which was launched on 19 November 2016, is the first satellite in the GOES‐R series. NOAA declared GOES‐16 to be operational on 17 December 2017. According to the Washington Volcanic Ash Advisory Center (VAAC), the 7 October Popocatépetl cloud extended to 30,000 ft (9.1 km) and the 4 November cloud reached 25,000 ft (7.6 km). In addition, two Sheveluch events were detected using 10‐min Himawari‐8 data. The advisories, issued by the Tokyo VAAC, indicate that the cloud detected on 10 October (23:40 UTC) reached 30,000 ft (9.1 km) and the cloud detected on 2 November (05:10 UTC) had a maximum height of 26,200 ft (8.0 km). Finally, four Ambae events were detected using the CGA method and 10‐min Himawari‐8 data: 13 October (07:00 UTC), 13 October (22:00 UTC), 22 October (3:10 UTC), and 29 October (04:00 UTC). The Wellington VAAC indicated that the 13 October Ambae clouds went up to 30,000 ft (9.1 km), the 22 October event generated a 12,000 ft (3.7 km) plume, and the 29 October cloud had a maximum height of 20,000 ft (6.1 km). The eight confirmed CGA alerts in the 6 October to 6 November 2017 timeframe captured the eruptive clouds in the very early stages of development—prior to exhibiting robust spectral signatures associated with the presence of ash and/or SO 2 . The z‐scores of the confirmed alerts ranged from 6.3 for the 4 November Popocatépetl case to 18.4 for the 10 October Sheveluch case. The typical latency of the alerts varies from 5 min (GOES‐16) to 26 min (Himawari‐8), depending on the sensor. Most of the latency is due to satellite data acquisition, not algorithm processing. Based on the Global Volcanism Program weekly reports and VAAC advisories, one volcanic cloud that reached 30,000 ft was missed by the CGA method in the 6 October to 6 November 2017 period. The 20 October 2017 (19:20 UTC) event at Tinakula (Solomon Islands; Global Volcanism Program, 2017) went undetected by VOLCAT because the cloud tracking failed the quality control test described in section 3.5. Of all of the events analyzed, including through near‐real‐time processing, this is the only known explosive event that was missed due to the quality control procedure. While some improvements are certainly possible, the CGA logic has been shown to be consistently effective at detecting volcanic clouds that reach the upper troposphere and beyond. Some lower level eruptions are also captured. Beyond the examples presented in this paper, many other confirmed alerts have been generated in near real time by VOLCAT over the last several years. A more comprehensive review and validation of the VOLCAT alerting capability is the subject of ongoing research. As stated earlier, the CGA‐based alerting tool does generate some false alerts. Figure 11 shows the number of false alerts per day during the 6 October to 6 November 2017 period as a function of satellite and VAAC region of responsibility (http://gis.icao.int/VOLCANO2014/). In this analysis, an alert is taken to be false if no volcanic ash advisory, consistent with the detected cloud, was issued. This is reasonable since volcanic ash advisories are based on several sources of information, including volcano observatories. While not always detected in the early stages (hence the motivation for an automated alerting capability), explosive events (hopefully) rarely go completely undetected. Further, our manual examination of the satellite imagery, associated with the false alerts, did not give reason to doubt the false alert classification. Figure 11 shows that the CGA method produces zero to eight false alerts per day, when the results from the GOES‐13, GOES‐15, MSG‐3, and Himawari‐8 satellites are aggregated. The combination of GOES‐13, GOES‐15, MSG‐3, and Himawari‐8 provides near‐global geographic coverage. In the median, three false alerts are generated each day. For additional context, the total number of cloud objects, tracked between consecutive satellite images, each day, is shown in Figure 11. The subset of those cloud objects, which exhibited behavior consistent with vertical growth, is also displayed in Figure 11. Approximately 1.3–1.6 million cloud objects are analyzed for vertical growth each day, with about 500,000 having a decrease (increase) in BT min (ε tot_max ) from one image to the next, consistent with vertical growth. Thus, the median false alert rate is quite small (approximately two false alerts for every million cloud objects analyzed or less than 1 out of every 100,000 objects exhibiting signs of vertical growth). The Himawari‐8 satellite produces the greatest number of false alerts because the data volume is much larger than the other satellites. When satellite‐specific object counts are accounted for, the performance (false alert rate) of each sensor is comparable. The number of false alerts per VAAC region is a function of the satellite coverage and the geographic extent of the VAAC region. The Tokyo VAAC region, which is imaged by Himawari‐8, has the greatest number of false alerts (total of 26). In the worst case, the 32‐day total of 26 alerts in the Tokyo VAAC region equates to less than one per day. Thus, the number of false alerts is manageable. The GOES‐16 satellite is not included in the false alert analysis since it was not yet operational during the selected analysis period. However, early indications are that the false alert rate is comparable to, or perhaps lower than, Himawari‐8. The GOES‐16 results will be reported in a future paper. Figure 11 Open in figure viewer PowerPoint The number of alerts from 6 October to 6 November 2017, not associated with known volcanic activity, is shown as a function of (top) satellite and (bottom) Volcanic Ash Advisory Center (VAAC) region. The total number of image‐to‐image cloud object tendencies that were analyzed (black sold line) and the subset of those tendencies exhibiting behavior that is consistent with vertical growth between two successive satellite images (gray solid line) are also shown on each panel relative to the right y axis, with units of millions. Each false alert in the 6 October to 6 November 2017 period was examined in detail. While it is occasionally difficult to completely rule out eruptive processes, as an influence on cloud development, from examination of the satellite data alone, human expert interrogation of satellite image sequences is usually effective for classifying CGA based alerts. Human experts can readily identify false growth due to tracking errors or measurement artifacts. In addition, human experts can often distinguish between volcanic and nonvolcanic cloud development by noting nuances associated with cloud origin and development (such as consistency with meteorological conditions) that are more difficult to characterize with computer algorithms. The manual analysis of false alerts indicated that about 60% of the false alerts can be attributed to cloud tracking errors, with the other 40% representative of actual changes in cloud object properties. Not surprisingly, the majority of tracking errors were associated with complex, fast‐moving, high‐level clouds. In addition, even when clouds are correctly tracked, the change in BT min and ε tot_max is sometimes due to advection over a colder background as opposed to actual vertical development. Fortunately, human experts can often identify the advection effect using the satellite imagery provided with the alerts. Future improvements will focus on reducing the number of false alerts as much as possible.

7 Conclusions A unique new geostationary satellite algorithm for automatically detecting emergent volcanic clouds, produced by explosive events, was described and demonstrated. The algorithm, referred to as the Cloud Growth Anomaly (CGA) method, identifies cloud objects that exhibit anomalous vertical growth between two successive geostationary satellite images, compared to meteorological clouds. While the CGA technique can be applied to any geostationary satellite sensor, as long as the image refresh time does not exceed 60 min, it is most effective when applied to latest generation of meteorological satellites such as Himawari‐8 and GOES‐16, which provide more frequent images with improved spatial resolution. The CGA method adds considerable value to the volcanic cloud detection problem because it does not rely on unique spectral signatures caused by differential extinction of radiation due to the presence of ash and/or SO 2 . Unique spectral signatures are often absent or difficult to conclusively identify in emergent volcanic clouds because of cloud opacity effects and/or the presence of hydrometeors (mainly ice). The CGA method was used to show that the very early development of volcanic clouds, produced by a large and diverse sampling of explosive events, is characterized by anomalous vertical growth, thereby allowing thresholds for automated volcanic cloud detection to be set. The automated detection capability is used within a volcanic eruption alerting system (VOLCAT), which has been running in near real time since 2013. The CGA‐based alerts, which are experimentally distributed to several Volcanic Ash Advisory Centers, have detected numerous eruptions, while resulting in less than one false alert per VAAC region, per day, on average. Volcanic cloud alerts are an important supplement to routine human expert analysis of satellite imagery, as rapidly increasing satellite data volumes do not allow for detailed manual analysis of every satellite image relevant to volcanic cloud detection. While the CGA method is not designed to detect all volcanic clouds, it can be combined with other detection techniques (e.g., Brenot et al., 2014; Pavolonis et al., 2015a, 2015b) to provide more complete, timely, detection of eruptive events, which is critical for aviation safety. Possible enhancements to the CGA algorithm include incorporation of localized trends in meteorological convection and satellite‐inferred surface temperature trends near volcanic vents (e.g., Wright et al., 2004). The CGA method has several potential uses beyond eruption alerting. Pouget et al. (2013) showed that time trends in certain spatial properties (e.g., cloud area) of volcanic umbrella clouds and vigorous volcanic plumes can be used to estimate mass eruption rate, which is a crucial eruption source parameter required by dispersion and transport models (Mastin et al., 2009). The CGA algorithm, in combination with other automated feature tracking methods, can be used to determine the cloud properties required to estimate mass eruption rate over time. While the cloud tracking technique employed with the CGA algorithm is not designed for birth‐to‐dissipation volcanic cloud tracking, it is important for the initial detection required for performing such long‐term automated tracking. Combining the CGA volcanic cloud detection capability with techniques better suited for middle and late life cycle tracking is the subject of ongoing research. In addition to determining eruption source parameters, automated volcanic cloud tracking is important for other modeling applications, such as supporting systematic verification and alternate methods of model initialization (e.g., Crawford et al., 2016). Finally, integrating satellite algorithms, like the CGA, with other measurement capabilities, including lightning detection, volcanic infrasound, and seismic monitoring (Fee & Matoza, 2013; Schneider et al., 2017; Van Eaton et al., 2016), is critical for improving volcanic cloud detection and characterization, as no single‐measurement source completely captures all aspects of volcanic cloud formation and evolution.

Acknowledgments The NOAA GOES‐R Risk Reduction Program funded this research (Federal Grant number: NA15NES4320001). The geostationary satellite data utilized in this study are available from the University of Wisconsin Space Science and Engineering Center (SSEC) Data Center (http://www.ssec.wisc.edu/datacenter/). The views, opinions, and findings contained in this report are those of the author(s) and should not be construed as an official National Oceanic and Atmospheric Administration or U.S. Government position, policy, or decision.

Appendix A: Cloud Object Tracking Search Box t = t 2 (hereafter referred to as the object of interest) is present in the t = t 1 image, a search box is defined in image coordinates. The starting and ending pixel coordinates of the search box (X 0 , X 1 , Y 0 , and Y 1 ) is a function of the difference between the two successive satellite images ( Δt ), the mean horizontal dimension of the satellite pixels ( ), the minimum and maximum x and y image coordinates of the object of interest (X min , X max , Y min , and Y max ), and an assumed wind speed ( |V| ). Thus, for a given pixel resolution, the size of the search box, in image coordinates, increases as a function of the time difference. The wind speed ( |V| ) is taken to be 100 m/s for objects with at least one cirrus‐like pixel (see Appendix x coordinate (scan line elements). The starting and ending y coordinate (scan lines) are computed in an analogous manner. In equation Δp is the number of pixels (rounded to the nearest integer), along the x or y direction, that a cloud feature will be displaced by the assumed wind speed. In equation N x is the maximum valid index. No assumptions are made regarding the wind direction, as the box accounts for movement in any direction. The search box is generally larger than needed and need not be precise since it is only used to improve the computational efficiency of the cloud object tracking algorithm. (A1) (A2) (A3) (A4) In order to efficiently determine if a given cloud object at time(hereafter referred to as) is present in theimage, a search box is defined in image coordinates. The starting and ending pixel coordinates of the search box (, and) is a function of the difference between the two successive satellite images (), the mean horizontal dimension of the satellite pixels (), the minimum and maximumandimage coordinates of the object of interest (, and), and an assumed wind speed (). Thus, for a given pixel resolution, the size of the search box, in image coordinates, increases as a function of the time difference. The wind speed () is taken to be 100 m/s for objects with at least one cirrus‐like pixel (see Appendix B ) and 50 m/s otherwise. Equations A1 through A4 illustrate how the starting and ending pixel coordinates of the search box are computed for thecoordinate (scan line elements). The starting and endingcoordinate (scan lines) are computed in an analogous manner. In equation A1 is the number of pixels (rounded to the nearest integer), along theordirection, that a cloud feature will be displaced by the assumed wind speed. In equation A4 is the maximum valid index. No assumptions are made regarding the wind direction, as the box accounts for movement in any direction. The search box is generally larger than needed and need not be precise since it is only used to improve the computational efficiency of the cloud object tracking algorithm.

Appendix B: Identification of Thin Cirrus Clouds Optically thin cloud regions, within cirrus shields, are difficult to track because they are often spatially expansive, fast‐moving, and/or transient. As such, a screening procedure is applied to limit the number of thin cirrus cloud elements that are tracked. The cirrus screening procedure improves the computational efficiency of the tracking routine and reduces the number of false cloud vertical growth signatures due to object tracking error. As described in Pavolonis (2010), Numerical Weather Prediction (NWP) model data and a fast radiative transfer model are used to simulate top of atmosphere infrared radiances under a variety of conditions. The NWP level where the emission from a simulated blackbody cloud and overlying clear atmosphere matches the observed radiance in a given channel can be readily determined. The highest tropospheric level, where the blackbody cloud assumption matches the observed radiance, is found for an infrared window and water vapor absorption channel. Every geostationary meteorological satellite has such spectral channels. Due to the absorption by water vapor, the atmospheric weighting function of a channel in the 6–7‐μm region will peak at much higher altitudes than a channel in the 10–12‐μm window region, where absorption by water vapor and other atmospheric gases is weaker. As such, for nonopaque clouds, the altitude at which the blackbody cloud assumption matches the observed radiance will be higher when computed using the water vapor absorption channel compared to an infrared window channel. Another way to capture the difference in altitude is through the NWP temperature associated with the simulated blackbody cloud level. When optically thin cirrus clouds are present, the difference in the temperature of the simulated infrared window and water vapor absorption channel blackbody clouds will be large, whereas the temperature difference will approach zero when an actual opaque cloud is actually present. This blackbody cloud differencing technique is analogous to simply taking the difference in the measured brightness temperature between an infrared window and water vapor absorption channel but is less sensitive to viewing angle and the channel spectral response functions, making application of a single threshold more effective. A heuristically determined threshold of 10 K is applied (e.g., >10 K implies that thin cirrus is present). Water vapor channels centered near 6.7–6.9 μm are utilized for all sensors, except MSG SEVIRI, where the water vapor channel utilized is centered near 6.2 μm (MSG SEVIRI does not have a channel 6.7–6.9‐μm channel). The infrared window channel utilized is centered near 10–11.5 μm. On sensors with two channels in this region, the channel with a center wavelength closer to 10 μm is utilized. In addition to the 10‐K threshold, the window channel cloud emissivity referenced to the top of the troposphere (ε tot ) must also be less than 0.75. All cloud objects where more than 50% of the pixels meet the thin cirrus criteria at t 2 (the second satellite image in the image pair used by the tracking algorithm) are flagged as thin cirrus and are not tracked.