Determining the height of deep volcanic eruptions over the tropical western Pacific with Himawari-8

2.1. Satellite

Satellite data used in this study come from the Japanese Meteorological Agency’s Himawari-8 satellite in geostationary orbit (Bessho et al. 2016). The primary instrument on Himawari is the Advanced Himawari Imager, a 16-channel instrument measuring radiance from the visible to the thermal infrared (IR). In the IR channels, the pixel resolution at the sub-satellite point is 2 km. For this work, channel 14 centred at 11.2 μm is the primary channel. Data are processed and archived by the Australian Bureau of Meteorology and accessed using Australia’s National Computing Infrastructure supercomputer (Australian Bureau of Meteorology 2021).

Because the stratosphere is isothermal or warms with height, identifying the height of volcanic plumes that extend into the stratosphere is difficult. When plumes are confined to the troposphere, the minimum BT can be identified with reasonable assurance that this temperature represents the coldest (and therefore highest) cloud top. When plumes are stratospheric, an examination of the horizontal structure of the volcanic cloud top is necessary to ascertain the height of the plume. Generally speaking, volcanic plumes that ascend into the stratosphere appear either as (1) persistent clouds (~1–12 h) that have a distinctive ‘stratospheric warm spot’ (SWS) signature or as (2) short-lived overshooting tops (OTs) that are notably colder than the surrounding cloud top. Conceptually, plumes that appear as an SWS are assumed to have stopped their rapid ascent, produced or transported enough ash, ice and cloud to produce a radiometric signal, and are in thermodynamic equilibrium with the ambient environment. Clouds that manifest as OTs are presumably regions of recent strong ascent, with largely undilute central updrafts. A parcel in this circumstance is negatively buoyant and is decelerating and beginning to mix with the ambient environment. As a result of these processes, OTs can evolve into SWSs on short time scales.

The first concept is illustrated in schematic form (Fig. 1), with colours representing BT contours as they might be visualised in satellite imagery. In the idealised tropospheric cloud (left), the coldest cloud portion of the cloud top is in the interior, and the gradient of the brightness temperature (∇BT) is axisymmetric about the coldest, highest portion of the cloud, with vectors directed outward towards the edges of the cloud feature. With the idealised stratospheric cloud (right), the contours in the centre of the cloud ‘reverse’ with warmer temperatures observed inside the coldest contour, creating a distinct signature, the SWS noted earlier. In many cases, only the SWS corresponds to the stratospheric cloud; most of the cloud remains in the troposphere as part of the umbrella. Corresponding to the SWS, ∇BT changes direction and locally points inwards towards the warmer cloud, creating divergence of the vectors in the ∇BT field.

Fig. 1. 

Schematic diagram illustrating behaviour of brightness temperature in tropospheric and stratospheric cloud tops. Colours represent infra-red brightness temperature contours, with red > yellow > green > blue, arrows depict the direction of the local gradient of the brightness temperature.

In many SWS cases, a stratospheric cloud can be readily identified via visual inspection of satellite imagery. The ability to do so depends on factors like the choice of enhancement used to visualise the imagery and the relative strength of the SWS, i.e. the BT contrast. In other events, the SWS is more subtle and difficult to identify. Following the discussion above, one useful metric is to evaluate the Laplacian of BT:

The calculation is performed at the scale of individual pixels using a standard centred finite difference approach. From a practical viewpoint, the pixel level calculations are noisy and can pick up on small fluctuations in BT; to alleviate this issue the results are smoothed over a 3 × 3 pixel area, putting the focus on larger-scale structures in the cloud top. Negative values indicate a reversal of the temperature gradient and possible SWSs as depicted Fig. 1. This calculation provides a useful guide to assist in the analysis; it does not exclusively and objectively identify SWSs without human intervention.

For the OT scenario, we adopt the approach described by Griffin et al. (2016). This approach calculates overshooting top height Z_OT:

where Z_U is the height of the umbrella (estimated from the umbrella BT, BT_U), BT_OT is the BT of the OT, and LR is the lapse rate, describing how the cloud cools above the umbrella/anvil. From a purely parcel theory perspective, this rate would be equivalent to the adiabatic LR; given the minute amounts of water vapour at the altitudes under consideration here, this is essentially the dry adiabatic LR = 9.8 K km⁻¹, modified by the physical process of entrainment. Owing to issues of sensor spatial resolution creating warm biases in the satellite imagery, Griffin et al. (2016) set LR = 7.3 K km⁻¹ for 1-km pixels from the Aqua satellite; lower spatial resolution onboard geostationary creates even larger biases; linear regression identified the magnitudes of these biases for 3- and 4-km resolution sensors (SEVIRI and GOES-13 respectively). Bedka et al. (2018) adopted this approach for the Australian region using the 4-km JMA MTSAT series of satellites, suggesting LR = 3.4–4.9 K km⁻¹. Using typical BTs for umbrella and OT observed here, the SEVIRI regression suggests a LR bias of ~1 K km⁻¹. The resolution of AHI is higher than SEVIRI, and the modelled bias may be an overestimate. For this work we assume LR = 6.5 K km⁻¹, equivalent to a bias of 0.8 K km⁻¹.

2.2. Radiosondes

To create of standardised reference profile for the UTLS in the region of interest, radiosonde data from Version 2 of Integrated Global Radiosonde Archive (IGRA2; Durre et al. 2018) are used. Radiosondes are widely spaced, but frequently launched (sub-daily in many cases) with an extensive historical record. Profiles of pressure, temperature, humidity and wind are reported on pre-defined ‘mandatory levels’ at fixed pressures and on sporadic ‘significant levels’, points where the wind, temperature or humidity observations deviate from the prevailing tendencies of the profile in question. Radiosondes reliably reach altitudes of 30–35 km, well into the stratosphere. Radiosondes may be subject to homogenisation issues related to instrumentation changes and reporting techniques (e.g. Jovanovic et al. 2017). These are not accounted for in this work.

Given the limitations of radiosondes, data from multiple stations over an extended time period are required to construct an adequate reference profile. A further constraint is that the area from which the stations are drawn should show a high degree of spatial homogeneity. This study focuses on a large area of the tropical western Pacific, approximately 95–180°E to 15°N–20°S (Fig. 2) with 17 long-term stations chosen for analysis (Table 1). This area broadly covers tropical portions of the VAAC Darwin, VAAC Tokyo and portions of VAAC Wellington and VAAC Washington. Data are analysed over the period of 2000–2021 to ensure that the reference profile reflects observed stratospheric temperature changes of 1–3 K since 1979 (Steiner et al. 2020; Scherllin-Pirscher et al. 2021). At a minimum, all selected stations have records that extend back to the early 1970s.

Fig. 2. 

Map of 17 long-term radiosonde stations used to create the UTLS reference profile.

The aim of the analysis is to create a high-resolution (0.5 km) reference profile that can be used to convert satellite-based IR BTs into reliable, consistent estimates for the height of volcanic plumes that reach the UTLS. We confine the analysis here to pressures below 200 hPa, ~12.5-km altitude, using ‘mandatory’ and ‘significant’ level data for temperature (T), pressure (P) and geopotential height (Z). The Z is generally only reported on mandatory levels in IGRA2 data and so it is calculated using the hypsometric equation, neglecting water vapour (virtual temperature T_V ≈ T at these levels) and the Eötvös effect (Ong and Roundy 2020). The reported 200-hPa value of Z is used as the baseline for the calculation. The values of Z, P and T for both the ‘lapse-rate tropopause’ (LRT), based on the standard World Meteorological Organization (1957) definition, and the ‘cold point tropopause’, the coldest point in the sounding, are also estimated for each sounding; the reported LRT ‘flags’ in the data are often not reliable or missing.

The processed sounding data are placed into two joint frequency distributions (JFDs), one each for Z–T and Z–P. These JFDs have bin sizes of ΔZ = 0.5 km, ΔT = 1.0 K and ΔP = 1 hPa and cover the ranges Z = 12–42 km, T = 178–278 K and P = 1–200 hPa. The JFDs are tabulated over months, with all years in the analysis period included in the appropriate month. Seasonal and annual means are determined by combining JFDs for individual months. Each of the 17 stations has its own set of JFDs. To generate the profiles, means are computed along height levels of the JFD to get the pressure and temperature profile at a given station; uncertainty is estimated by determining the range from the 16th to the 84th percentile at a height level, equivalent to a 1 standard deviation (s.d.) span where the data are normally distributed. All-station means for the larger area are estimated by averaging profiles generated at the individual stations.

2.3. Supplemental data

Additional data from several sources are used in the discussion of the environments of the 19 December 2021 eruption of HTHH. First, profiles of European Centre for Medium-range Weather Forecasts ERA5 reanalysis (Hersbach et al. 2020) T and Z interpolated to pressure levels with the standard 0.25° horizontal resolution within a 6 × 6° box centred on the volcano are used. Second, Global Navigation Satellite System Radio Occultation measurements made from the University Corporation for Atmospheric Research (UCAR) Constellation Observing System for Meteorology Ionosphere and Climate (COSMIC-2) are used (UCAR COSMIC Program 2019). Specifically, the COSMIC-2 ‘wetPrF’ 1D-var retrieval data within 1000 km and ±6 h of the eruption time are used.

3.1. Joint frequency distributions

The range of data quality and sampling is seen by comparing the Z–T JFDs for the radiosonde stations at Yap and Medan (Fig. 3). Yap is an excellent sample, whereas that at Medan is weaker, but still adequate for the task. From the number of observations in each height bin it is obvious that the height bins that contain mandatory levels are better sampled, with 2–3 times more observations than the others. The better sample at Yap has a larger number of observations extending up to 35-km altitude and beyond, well above the 10-hPa mandatory level, whereas at Medan the total number of observations drops off above ~25-km altitude.

Fig. 3. 

Annual geopotential height–temperature (Z–T) joint frequency distributions (JFDs) at Yap (left) and Medan (right). Shaded contours represent the frequency observed in each bin, with contours at [1e−5, 1e−3, 1e−2, 3e−2, 5e−2, 7e−2, 1e−1, 2e−1]. Thick yellow line represents the mean profile, with horizontal bars being the approximate 1-s.d. spread. Thick black contours represent the JFD for LRT height, with contours at [1e−4, 1e−3, 1e−2, 3e−2, 5e−2]. Red plot to right of each panel is the total number of samples in each height bin.

In both examples, the data are well clustered around the mean, with the middle two-thirds of the data confined to within 3–4 K of the mean, but with large tails. The distributions are narrower in the troposphere, but wider in the stratosphere where the gravity waves dominate. The largest variability is near the tropopause at 16–18 km, where temperature range is 180–205 K. Tropopause heights and temperatures vary similarly; at the more poorly sampled stations, the tropopause is more erratic with more frequent outliers, but still generally consistent. Overall, although there is variability in the sampling of the JFDs, the number of observations are generally adequate below heights of 30–35 km and the reliability and spread of the data are adequate for the task.

3.2. Annual and seasonal means

Station profiles from each of the 17 stations are computed and all are averaged to all-station mean pressure and temperature profiles representative of the entire region. The annual all-station mean Z–P profile (Fig. 4a) shows remarkable consistency with height across stations. The logarithmic variation of pressure with height is evident. Plausible, consistent pressure estimates are available up to 38-km altitude. The annual all-station mean Z–T profile (Fig. 4b) also shows strong consistency across stations, with a reliable temperature profile available up to ~35 km. The largest temperature spread is identified near the 17 km, with a second spread in the data between 20 and 25 km. Although there is more variability between stations compared to the pressure profile, there is a remarkable spatial consistency in the data. The consistency of the results between stations across the network builds confidence that a suitable reference profile is being captured in the analysis.

Fig. 4. 

Annual all-station mean (a) Z–P and (b) Z–T profiles (black and yellow dashed line) with annual means from individual stations (thin coloured lines).

The complex structure of the TTL is well-captured by this analysis. The main portion of the TTL here is identified within ~16–18 km, a region where vertical temperature variation is minimal as the LR transitions from cooling-with-height to warming-with-height. The analysis also captures the ‘tropopause inversion layer’ in the lower stratosphere (e.g. Bell and Geller 2008), a region of enhanced static stability above the TTL, despite not being in ‘tropopause-relative coordinates’. Overall, the general temperature structure approximately matches in absolute value with the profile (up to 25 km) presented in fig. 1 of Scherllin-Pirscher et al. (2021).

Profiles based on the traditional meteorological seasons (e.g. December–January–February, DJF; March–April–May, MAM) reveal the magnitude and direction of the annual cycle (Fig. 5), with the upper and lower stratosphere behaving differently across the year. The TTL (and the stratosphere below ~22 km) is at its highest and coldest during DJF, and at its warmest and lowest during JJA. Conditions are close to the annual average during the equinoctial seasons. Above 26 km, the temperature profile is warmest during MAM and coldest during JJA and DJF. There is little seasonal variation in the upper troposphere below 16 km.

Fig. 5. 

Z–T profiles for the annual (ANN) and seasonal (DJF, MAM, JJA and SON) all-station means.

Both the direction and magnitude of these seasonal variations agrees well with the monthly variations of stratospheric temperature from ERA-40 reanalysis of Fueglistaler et al. (2009). That study identified a peak-to-peak change in temperature of ~8 K at the 80-hPa level; in this work, the maximum T difference is ~7 K at 17.75 km, ~81 hPa. In the upper stratosphere, anomalies in T from the annual mean are ~1.5 K in MAM and cooler in other seasons, again in good agreement with Fueglistaler et al. (2009). The seasonal variation in the Brünt–Vaisala frequency and associated structure of the tropopause inversion layer (not shown) agrees well with Scherllin-Pirscher et al. (2021).

The seasonal uncertainty profile estimates for the means are useful for understanding the vertical structure of the uncertainty (Fig. 6). The numbers represent the station averaged 16–84th percentile range divided by 2 to approximate the s.d. In the upper troposphere, below 16.5 km, uncertainty is typically 2–3 K. The TTL, ~16.5–20 km, is where the largest uncertainty is found. Here it ranges within 3–5 K, with the largest uncertainty noticeably occurring in DJF. In the stratosphere (20–35 km), the uncertainty is 3–3.5 K at most heights.

Fig. 6. 

Seasonal profiles of temperature uncertainty.

3.3. Seasonal reference profiles

The consistency of the results and the strong agreement with previous studies indicate that the technique to construct the profiles is appropriate. The structure and variability of the UTLS annual cycle over the tropical western Pacific is adequately captured in this analysis; other sources of the multi-year variability driven by the QBO and ENSO and the short-term gravity wave fluctuations are inherent in the uncertainty estimates of Fig. 6. With that in mind, Table 2 presents standard reference profiles for the region. The left column is the height, with the numbers under the middle columns being the reference temperature at that height. The right column is an estimate of the uncertainty of the height estimate, derived from the horizontal uncertainty (Fig. 6) and the vertical structure of the profile by considering the change in height resulting by adding or subtracting the uncertainty value to or from the reference temperature in that layer. For all seasons, the profiles are extrapolated to a nominal stratopause by assuming a constant LR of +2 K km⁻¹. This is equal to the mean LR for 25–35-km altitude in the profiles, approximately equal to the LR in the US Standard Atmosphere (Minzner 1977), but with different starting values.

4.1. 20:40 hours UTC image

Fig. 7 shows an image of the eruption from the Himawari-8 scan beginning at 20:40 hours UTC on 19 December 2021, the first satellite scene where the eruption is visible. The actual sampling time at the latitude of the volcano is ~7.5 min after the nominal image start time. The image combines a Rayleigh-corrected ‘True Colour’ RGB (e.g. Broomhall et al. 2019) with selected contours of AHI channel 14 (11.2 μm) IR brightness temperature and negative ∇²BT, indicating a reversal of the gradient. The scene is very complex and tells a significantly different story to the official account.

Fig. 7. 

Rayleigh corrected 1000-m ‘True Colour’ imagery of the Hunga Tonga-Hunga Ha’apai eruption from the Himawari-8 scan beginning on 20:40 hours UTC on 19 December 2021. Black contours show 2000-m channel 14 (11.2 μm) BT between 200 and 215 K at 5 K intervals. Green contours are ∇²BT at values of −0.1 (solid) and −1.5 (dashed) indicating a reversal of the cloud top temperature gradient. Latitude–longitude grid lines are spaced every 0.1°. The nominal location of the volcano is noted with the (red) triangle. Annotations are discussed in the text.

The cloud top portion annotated with a ‘1’ is most reminiscent of a traditional anvil or umbrella. The coldest BTs here are 201–202 K; the reference profiles (Table 2) indicate that this region lies below the tropopause, 15 ± 0.5-km altitude. This value agrees with the analysis of the umbrella height for this event by Gupta et al. (2022). Near ‘2’, a significant SWS is noted by the negative ∇²BT and is quite obvious in the IR contours. Note that the smoothed ∇²BT does not fully capture the extent of the SWS; this is a function of the smoothing as the strong outward pointing gradient to the west nearer the cloud edge ‘corrupts’ the score. Maximum IR temperatures in the SWS are 218.8 K. Also note the dark spot near ‘2’, which is at least partially due to a shadow cast by the ‘brighter spot’ at ‘3’. The BT at this cloud feature is 208–209 K, with Table 2 suggesting a height of 22 km. Bright cloud tops noted at ‘4’ and ‘5’ vary in BT between 206 and 214 K, indicating that the plume top height is 21–24 km. Also note the weak OT near ‘6’; its minimum BT is 197.9 K, with the OT equation indicating that it is ~500 m higher than the umbrella top at ‘1’, ~15.5 km, in general agreement with what would be predicted using Table 2.

In this scene, there is strong evidence that the eruption cloud top height reached at least 24 km in the first 10 min of this eruption. The warmest cloud top feature is at ‘2’ and suggests the height at least briefly reached 26 km from Table 2. Uncertainty in these heights estimates is ±2 km. The shadow or dark cloud does not indicate an active, growing plume at that location. Careful examination of the image does suggest a subtle shift in colour to the surrounding background. One hypothesis is that the eruption grew to that height, but collapsed quickly due to a loss of buoyancy, leaving only a transient cloud that has detrained from the plume. A second possibility is this is a pileus cloud generated above the rapidly rising plume. Given the dryness of the stratosphere in general, the first possibility is more probable. In either case, these clouds must be optically thick enough to generate the warm temperatures observed by the satellite.

4.2. 20:50 hours UTC image

Fig. 8 is an image in the same format as Fig. 7, but at 20:50 hours UTC; as before, the scan time of the volcanic cloud is 7.5 min past the nominal image time. The volcanic cloud has increased in area by a factor of ~2.5 from the previous image; the area covered by the SWS has similarly grown, and now contains three areas of stronger ∇²BT. As before, the identification of the obvious warm spot is not perfect as the smoothing of regions near the cloud edge dilutes the signal. The umbrella (‘1’) remains at ~201–202 K as it moves off to the north and north-east. Table 2 suggests that this is near 15-km altitude. The regions marked ‘2’ are the locations of the high cloud tops annotated as ‘4’, ‘5’ and ‘6’ in Fig. 7. The rightmost one aligns with the incipient OT noted previously at 15.5-km altitude (‘6’); it now has an IR BT of 206 K and ~21-km altitude based on Table 2. This illustrates the transition of an OT into an SWS over a 10-min period. The other regions marked ‘2’ (previously ‘4’ and ‘5’) have BTs that are slightly cooler than in the previous image, suggesting the previously identified clouds are not growing vertically or have begun to collapse.

Fig. 8. 

As in Fig. 7, except at 20:50 hours UTC 19 December 2021. Note that area covered by the image is approximately twice that of Fig. 7.

The region marked ‘3’ is the warmest part of the cloud top with several pixels with IR ~213 K, ~23 km based on Table 2. Close examination of the True Colour image indicates a different colour and texture of the cloud here. To the east of this feature, a significant OT cloud is apparent; a significant shadow is visible but is blocked by the higher cloud of the main plume. The minimum BT of 189.4 K indicates this extends 1.8 km over a 201 K anvil, ~16.8 km in height.

4.3. Later imagery

Summarising later imagery (not shown), an SWS continues to be maintained and replenished until at least 21:50 hours UTC in the vicinity of the volcano. Several further OTs are identified, and a maximum warm spot temperature of 219.8 K is observed at 21:30 hours UTC, placing the top above the 26 km level. After 21:50 hours UTC, an SWS is detectable, but decaying and drifting toward the south–south-east until 22:30 hours UTC. This analysis suggests that the eruption was well into the stratosphere for 1–2 h. After that time, cloud top BTs gradually increase and the cloud becomes optically thin. Simultaneously, the umbrella in the troposphere continues to spread to the north and north-east and reaches widespread minimum BTs of 195–198 K, placing it within ~15.5–16.5 km. The eruption continues although reaching much lower altitudes; BTs in the active plume do not go below 200 K from 01:00 hours UTC on 20 December 2021. The apparent eruption on IR ceases at 05:20 hours UTC, although the remnants of the volcanic cloud remain detectable beyond 10:00 hours UTC.

4.4. Choice of temperature profile

Fig. 9 shows temperature profiles from a variety of sources on 19 and 20 December. Aside from the climatological profile, the only information on this plot that was operationally available are the three radiosonde profiles from Nadi, Pago Pago and Funafuti, 800–1400 km distant from the eruption at the nominal time of 00:00 hours UTC on 19 December, nearly 1 day before the eruption (solid, coloured lines). Considerable differences are seen in the profiles, illustrating the difficulty in choosing a representative profile. Individual profiles, including those from ERA5 at the nominal eruption time, are noisy and show large variations both spatially and temporally; the influence of gravity waves is apparent in the vertical structure. The largest differences between data sources are within the TTL. Overall, the DJF climatological profile provides the best match to the detailed radiosonde observations at those heights. Some individual COSMIC-2 profiles capture the extremely cold temperatures, whereas the ERA5 pressure-level data poorly resolves the complex structure. Above the TTL, between 18 and 26 km, the ERA5 and COSMIC-2 mean profiles are up to ~4 K warmer than the climatological profile; individual profiles, including the radiosondes, vary considerably, being both warmer and cooler than the climatological profile. On average, height estimates are 1–2 km lower using data from the day, although those estimates both above and below that of the climatological profile all lie within the estimated uncertainty of the climatological profile.

Fig. 9. 

UTLS temperature profiles from different sources on 19 and 20 December 2021. Shown are the DJF climatological profile (red dash), the three closest radiosonde sites (green, orange and pink) on 19 (solid) and 20 December (dash–dot) at 00:00 hours UTC. Also shown are 16 COSMIC-2 profiles (cyan) and their mean (blue dash) and 625 profiles from the ERA5 reanalysis (grey ‘shading’) and their mean (black dash). The map inset shows the locations of Hunga Tonga-Hunga Ha’apai (red diamond), the radiosonde sites (squares), COSMIC-2 profiles (blue crosses) and the ERA5 box (grey).

This analysis shows that the climatological profiles are accurate enough for operational usage, particularly accounting for the largely unknown uncertainties associated with individual observations. For post-event analysis, heights of volcanic plumes may require refinement using data not available in real time; appropriate averaging should be used in these cases to remove the short-term effects of gravity waves in individual profiles, particularly in the stratosphere.