Stereo Plume Height and Motion Retrievals for the Record-Setting Hunga Tonga-Hunga Ha'apai Eruption of 15 January 2022

Stereo methods using GOES-17 and Himawari-8 applied to the Hunga Tonga-Hunga Ha'apai volcanic plume on 15 January 2022 show overshooting tops reaching 50–55 km altitude, a record in the satellite era. Plume height is important to understand dispersal and transport in the stratosphere and climate impacts. Stereo methods, using geostationary satellite pairs, offer the ability to accurately capture the evolution of plume top morphology quasi-continuously over long periods. Manual photogrammetry estimates plume height during the most dynamic early phase of the eruption and a fully automated algorithm retrieves both plume height and advection every 10 min during a more frequently sampled and stable phase beginning 3 hr after the eruption. Stereo heights with Global Satellite System Radio showing that much the


Materials and Methods
We rely primarily on automated stereo methods previously developed at NASA for tracking wind tracers from multiple satellites ("stereo winds"). These tools have already been applied to stereo observations from Himawari-8 and MODIS of the 2019 eruption of Raikoke Horvath, Girina, et al., 2021) and a similar approach has been applied to the study of an eruption of Mt. Etna in 2013 using a pair of Meteosat satellites (Merucci et al., 2016). Our stereo-winds code (Carr et al., 2020) uses a triplet of images from one geostationary satellite and imagery from another satellite after it has been remapped into the same geometry as the former. Tracers in the geostationary satellite imagery are tracked in time and between satellites by pattern matching. Our code offers the ability to retrieve both the height and the motion of a volcanic plume nearly continuously until the eruption dies down and the plume becomes too tenuous to track effectively.
We use geostationary satellite imagery from GOES-17 nominally stationed at 137.2°W and Himawari-8 nominally stationed at 140.7°E. Both satellites have nearly identical 16-channel imagers and capture the Earth's Full Disk (FD) every 10 min. Both imagers have similar capabilities to acquire smaller sectors with shorter refresh times. In particular, GOES can acquire an approximately 1,000 × 1,000 km 2 Mesoscale (MESO) scene every 1 min. NOAA centered a MESO scene over the eruption at 07:05Z to begin covering the volcano every minute. Before 07:05Z, the volcano was only covered in the FD scenes repeating every 10 min. Stereo-winds feature tracking works better when the time between scene repetitions is shorter. At the beginning of this eruption, as the radius of the plume top expanded at a mean radial speed of ∼50 m s −1 from ∼20 km at 04:10Z to ∼55 km at 04:20Z and to ∼80 km at 04:30Z, features that otherwise could be tracked morphed quickly in the billowing and rapidly expanding plume. This made automated feature tracking more difficult using the 10-min GOES-17 FD scenes. Therefore, we focused our use of the stereo-winds code on GOES-17 MESO scenes paired with a single near-simultaneous Himawari FD. For example, the Himawari FD beginning at 07:00Z, starts above the North Pole and progresses toward the south, reaching the volcano around 07:06Z; therefore, we pair GOES MESO scenes at 07:05Z, 07:06Z, and 07:07Z with the single 07:00Z Himawari FD scene. This pattern repeats every 10 min. The sun has already set at the volcanic site by 07:00Z, so we use the long-wave Band 14 (B14) with a center wavelength at 11.2 μm and spatial resolution of 2 km for both satellites.
Simultaneous stereo pairs are not required by the stereo-winds code as pixel acquisition times are used in the retrieval model. In our case, the 07:06Z GOES-17 MESO and the 07:00Z Himawari FD differ from exact simultaneity by either ∼10, ∼40, or ∼70 s across the mesoscale domain depending on which of the three 30-s swaths of the Himawari FD that cover the mesoscale domain contain the pixel. The apparent displacement of a feature between images represents a combination of motion and geometric parallax that the model separates knowing the pixel times. This is a key feature that enables diverse combinations of independently operated satellites to be used.
During the most dynamic early phase of the eruption before 07:05Z, however, we relied on manual photogrammetry Prata & Grant, 2001). The height of a feature visually matched in two images was determined from the ellipsoid-projected parallax between the match points and the corresponding satellite view/ azimuth angles. Feature matching was performed between the highest resolution (0.5 km) GOES and Himawari red band (0.64 μm) FD pairs with a typical non-simultaneity of ∼30 s in each 10-min slot. The method can also be applied to apparent shadows observed in a single image, in which case the second satellite view is effectively replaced by the solar-projected location of the feature, that is, the shadow terminus. Note that the method can estimate the height but not the motion of a feature.
Our stereo heights are confirmed by GNSS-RO profiles. The GNSS-RO technique has been used operationally in numerical weather prediction to provide temperature information in the lower stratosphere and upper troposphere through the bending induced by atmospheric refractivity (Healy et al., 2005). BA is the critical parameter that allows profiling atmospheric variability; it increases exponentially as the surface approaches. To detect eruption effects on BA, we compared the perturbed BA profiles against the background, using near real-time data from the Constellation Observing System for Meteorology, Ionosphere and Climate-2 (COSMIC-2), and Spire under NOAA's Commercial Weather Data Pilot (CWDP) program.  Figure 1 shows the maximum height estimates found at 04:30Z. The plume resembled a stack of pancakes with several identifiable layers, capped by a large central dome (see the contrast-enhanced Supporting Information S1). Heights were calculated for 30 plume retrievals and one low-level cloud at 2 km altitude serving as a reference. The stereo retrievals indicate layers at 20-21 km, 26-27 km, and the massive topmost ring of ∼80 km radius at 38-40 km. There are significant km-scale variations within this rugged topmost ring, demonstrating the sheer dynamism of the explosive rise. Most remarkably, the central dome reached an altitude of ∼55 km. The retrievals show a systematic increase in height as sampling ascends from the base to the dome's peak. An independent estimate of peak height from applying the side view method of  to the corresponding GEO-KOMPSAT-2A image yields ∼57 km, as described in Supporting Information S1. The GOES-17 view was favorable for applying the shadow method and is unaffected by motion effects. The long white arrow marks the apparent shadow of a local peak in the topmost ring, cast on the ocean surface. For a solar zenith/azimuth angle of 68°/−108°, the 145.7 km shadow length corresponds to a plume edge height of ∼40 km, in good agreement with the nearby GOES-17-Himawari-8 stereo estimates. The shorter white arrow indicates the shadow of a dome feature cast on the ring. This 61.7 km long shadow corresponds to an ∼18 km height differential relative to the ∼38 km altitude of the ring feature at the shadow terminus, consistent with a peak dome altitude of ∼55 km within the expected tolerance of ±2 km. Shadow-based plume height estimates for the 04:50-05:00Z time slots are also in Supporting Information S1.

Automated Stereo-Winds Method
The stereo-winds method was applied after 07:05Z when NOAA placed the GOES-17 MESO1 scene on the volcano. The 1-min cadence of the MESO scenes is nearly ideal as it allows enough time for motion to be observed but leaves feature shapes sufficiently invariant to be accurately tracked with subpixel precision. Figure 2a shows the jointly retrieved stereo heights, and advection vectors for 08:06Z (uses GOES-17 MESO scenes at 08:05Z, 08:06Z, and 08:07Z) paired with the Himawari 08:00Z FD. Smaller tracking templates allow for a higher resolution characterization of the plume geometry and dynamics, but this must be balanced against tracking accuracy. We selected 16 × 16 pixel templates (32 km feature spatial resolution) and oversampled the scene 4:1 (8 km sampling). A lofted stratospheric plume mostly above 30 km is seen over a lower cold umbrella at ∼18 km near the tropopause. The expansion and westward drift of the plume is evident from the advection vectors. Accompanying temperature assignments have been made according to the Effective Blackbody Temperature (EBBT) of the tracking templates taken from the 08:06Z MESO scene. Templates lower than 22 km (cold umbrella and troposphere) are presumed to be cold targets over warmer backgrounds and are assigned the mean temperature of the 20% coldest pixels. Templates 22 km and above are expected to be warmer than the cold umbrella and are assigned the mean temperature of the 20% warmest pixels because of the negative stratospheric lapse rate. However, where the plume becomes optically thin (effectively semi-transparent) near its edges, heat from the troposphere below makes the templates appear warmer than expected in the stratosphere up to 35 km (>240 K). In such cases, a colder pattern is being tracked over a warm background, and the 20% coldest pixels are averaged to be more representative of the plume temperature. This approach to temperature assignment is similar in concept to that proposed by Borde et al. (2014) where the assignment is being made based on the pixels that are most important for feature tracking.
ERA5 and GNSS-RO temperature profiles are provided in Figure 2b for context, which may not be entirely accurate in the presence of volcanic ash; however, the EBBTs of both the plume and tropospheric clouds broadly follow and cluster near both profiles. Deviations from the profiles are indicative of challenges with temperature-based height assignment, which in general is complicated by inversions, semitransparency, undercooling, temperature assignment uncertainty, and targets with emissivity less than one. Additional confidence in the stereo retrievals is gained by noting that the field of marine stratocumulus (Sc) clouds in the SW quadrant of the mesoscale domain has stereo heights between 500 m and 1.5 km. These are within expectations for such clouds and match well with the NOAA temperature-based cloud top height product for this scene.
The 08:06Z scene was selected for Figure 2 as an interesting case because it clearly shows an eruption of volcanic ash overshooting the plume top centered over the marked volcanic site (Global Volcanism Program, 2013). All stereo-winds product files are found in the Supporting Information S1, along with animations of the IR imagery, height-coded advection fields, and height versus temperature profiles for the complete unthinned set (Supporting Information S1). At 08:41Z, an overshooting plume core reached the upper stratosphere with an EBBT of ∼173K, significantly colder than the surrounding stratospheric plume. This cold bubble has an EBBT ∼15K colder than the lower umbrella near the tropopause, which can be explained by adiabatic cooling during plume convection (Woods & Self, 1992). Atmospheric gravity waves are also evident in the animations.
The evolution of the mesoscale plume height statistics from 07:06Z to 13:46Z for 15 January 2022 is detailed in Figure 3. This density plot, two-dimensional bivariate tiled histogram, presents the distribution of retrieved heights above 14 km binned over the mesoscale domain every 200 m and 10 min for nearly 7 hours. Optically  Panel (b) shows the assigned temperatures for each retrieval and the associated advection speed. Temperature profiles from ERA5 at 08:00Z (interpolated to the volcanic site) and COSMIC-2 RO at 07:11Z (23.34°S, 178.56°W) have been added for context. thick parts of the plume at a higher altitude can block the view of elements underneath, which explains the lower number density in the cold umbrella when the upper-level plume has a higher number density. If the upper volcanic plume were a mix of optically thin and thick layers, the height retrieval count profile would indicate an approximate upper-plume thickness of 4-5 km. However, as noted, optically thick parts will obscure material below. From 07:06Z to 10:06Z, the upper volcanic plume obscures the lower cold umbrella, decreasing its retrieval numbers, but as the upper volcanic plume moves out of the mesoscale domain, retrievals increase for the slower evolving cold umbrella. Subsequent activity is seen in the volcanic plume retrieval peaks between 07:00Z and 09:00Z that do not surpass the height of the initial eruption, for which earlier plume material has advected out of the mesoscale domain. This interpretation is supported by the animation of the stereo retrievals provided in the Supporting Information S1. Overall, the plume ceiling within the mesoscale domain is observed to fall nearly linearly at a rate of ∼1.7 km/hr.

Radio Occultation Results
To detect volcano-induced BA perturbations on January 15, we first derive a BA reference profile for the latitude band (15°S to 25°S) using all BA profiles from January 14 and 16. Then, we take the difference between individual profiles from January 15 and the reference BA to find anomalies near the Tonga eruption region. The BA perturbation is measured in percentage from the difference divided by the reference BA. Typically, a BA profile varies within 10% about the reference at heights between 20 and 40 km.
As shown in Figure 4, two strong BA perturbations are observed around 05:17:40Z (from COSMIC-2) and 05:12:42Z (from Spire), each about 150 km from the volcano, indicating an ∼50% deviation from the reference. The top heights of the BA perturbations are approximately 41 and 38 km, respectively for these cases. The GNSS-RO technique has a relatively good (<2 km) vertical but poor (∼300 km) horizontal resolution. Thus, the bending is sensitive to the refractivity change induced by a large volcanic plume like that of the Tonga eruption. The top of the BA perturbation indicates  the first height where the RO line-of-sight encounters the volcanic plume. Starting with the average of the two RO measurements, 39.5 km at 05:15Z, and extrapolating forward with the fall rate of 1.7 km/hr (identified in Figure 3) until 07:06Z predicts a maximum plume top of 36.4 km. This agrees with the highest retrievals in the 07:06Z histogram to the expected vertical resolution of the RO observations.

Discussion
Manual feature matching between the two satellite images is relatively straightforward; however, the horizontal and vertical motion of plume features introduces a height error due to the ∼30 s non-simultaneity between the views. For HT-HH, the GOES-17 and Himawari-8 view zenith/azimuth angles are around 50°/66° and 55°/−71°, respectively. This stereo geometry leads to a parallax almost precisely in the east-west direction, and a height error of ∼0.4 km for every 1 km parallax error (see Equation 8 in -the corresponding height error sensitivity for the shadow geometry is ∼0.3 km per 1 km parallax. In the early phases of the eruption, the velocity of the fastest moving plume element is ∼50 m s −1 in both the horizontal and vertical directions. In 30 s, this results in a ∼1.5 km horizontal movement and a ∼1.5 km vertical rise, the latter of which translates to a ∼3.7 km parallax increase. In the worst case, when these motion parallax errors add up, the total parallax error is ∼5.2 km, amounting to a ∼2 km height overestimation. Vertical velocity can be negative (but smaller in magnitude) when overshooting tops collapse. There can also be partial cancellation between the motion parallax errors, depending on the direction of the horizontal velocity. Therefore, a ±2 km height uncertainty is expected for the most dynamic plume features in Figure 1 (although overestimation is more likely than underestimation).
The automated height-resolved plume motion vectors are available every 10 min near real-time and thus have the potential to improve operational ash forecasts when assimilated into a dispersal model. The uncertainty in volcanic ash transport predictions arises from uncertainty in both the eruption source parameters (such as eruption column height) and the model wind field (Stefanescu et al., 2014). Recent work demonstrated a significant improvement of ensemble-based forecast skill when geostationary satellite observations of SO 2 mass loading were incorporated into the FALL3D atmospheric dispersal model Mingari et al., 2022). In a similar vein, our stereo height and motion retrievals could be used as the reference around which perturbations are calculated during the ensemble generation. This can be particularly useful for eruptions with a complex time-varying emission source term or when the observed atmospheric flow includes scales of motion unresolved by the model.

Conclusions
We used GOES-17 and Himawari-8 geostationary satellite imagery to study the 15 January 2022 eruption of the HT-HH volcano shortly after 04:00Z. Manually measured parallaxes between near-simultaneous visible images from each satellite and the lengths of shadows cast by the plume provided information about plume structure in the earliest phases of the eruption. We applied a stereo-winds tracking method starting at 07:06Z to automatically generate a high-density representation of plume advection jointly retrieved with plume height. This powerful method provided a detailed mesoscale description of advection and vertical structure for nearly 7 hr every 10 min.
Our results show that the HT-HH volcano ejected material as high as 55 km, with a substantial amount of material injected above 30 km. Such peak altitudes exceed the maximum column altitudes of 40 km reported for the 1991 Mount Pinatubo eruption by Holasek et al. (1996). As the eruption progressed and later died down, advection carried much of the injected stratospheric plume outside the mesoscale domain studied. The plume remaining in the mesoscale domain would have thinned or collapsed to render it untrackable in the IR imagery. Our analysis of GNSS-RO BAs confirms that material was lofted between 30 and 40 km approximately 1 hour after the initial eruption.