Modeling Volcanic Ash Dispersion from the Hunga Tonga-Hunga Ha'apai Eruption Using WRF-Chem and Meteorological FASDAS Data Assimilation

1. Introduction

Volcanic eruptions, because of the large amounts of gases and ash particles emitted into the atmosphere, exert a significant impact on the climate and the environment (Hirtl et al., 2019). Hunga Tonga Hunga Ha'apai is a submarine volcano located in the Pacific Ocean, about 160 kilometers northwest of the capital of Tonga, Nuku'alofa. The volcano erupted in December 2014 and again in March 2015, and it has been erupting periodically since then. The 2014-2015 eruption produced a new island, which has been called "Home Island" by the people of Tonga. The island has a size of about 5.5 square kilometers and is composed of volcanic ash and lava flows. The volcanic activity at Hunga Tonga Hunga Ha'apai has been monitored by the Tonga Met Service, the Pacific Tsunami Warning Center, and other agencies.

There have been no reports of significant impacts on people or infrastructure from the volcanic activity. However, the massive blast that rocked the island on January 15th 2022 was likely the largest one that caught the eye. The eruption was so violent that caused a tsunami that spread across the Pacific Ocean and whose shock wave circled the earth several times. Following the underwater volcanic explosion, dramatic satellite images showed the long, rumbling eruption of the Hunga Tonga-Hunga Ha'apai volcano spew smoke and ash in the air, with detected more than 190.000 lightning events during a 21-hour period between 14 and 15 January as recorded by the Global Lightning Detection Network (GLD360) ground-based network. In addition, significant ashfall was reported on populated islands of Tonga, which was the principal cause of fresh water supplies contamination and difficulty breathing from the ash in the air. It also caused flights to be cancelled and prevented marine transportation from and towards the Republic of Tonga.

The effect of the Hunga-Tonga Ha'apai volcanic eruption was clear on the state of the atmosphere in the troposphere layer, where many Surface Weather Observation Stations recorded a pressure disturbance at the moment of pressure wave pulses arrival (Yuen et al., 2022), which are known as Lamb waves (Amores et al., 2022). In recent study which conducted to more than 123 stations around the world, the propagation speed of the first Lamb wave was estimated at about 309±1 ${m.s}^{-1}$(Díaz & Rigby, 2022). For instance, after 14 hours pressure pulses were observed in more than 40 weather stations in the British and Irish Isles (Burt, 2022; Harrison, 2022).

On the other hand, although this noticeable effect is on the lowest layer of troposphere, the impact of this eruption on the upper atmosphere layers was greater. In this regard, Numerous studies indicated that the volcanic aerosols resulting from the eruption reached high altitudes in the stratosphere, up to 30 and 40 km (D’Arcangelo et al., 2022; Amores et al., 2022). (Kloss et al., 2022) indicated that the ash cloud eruption reached an altitude between 50 and 55 km. Generally, most of the studies currently conducted have focused on aerosols from equatorial volcanoes at altitudes between 20 and 30 km, due to aerosols are concentrated in the stratosphere between 20 and 25 km (Kremser et al., 2016). The aerosol remains for several years in the stratosphere (Robock, 2000).

The eruption occurred under the ocean and ejected a large amount of water vaper ($H_{2}O$) in the troposphere (Legras et al., 2022) and relatively small amount of ${SO}_2$ (Xu et al., 2022). Taha et al., 2022 estimated that the stratospheric aerosol optical depth from the eruption was comparable to that of Pinatubo eruption (1991), whilst the measurements showed that the water vapor continued to rise slowly after the explosion, and the aerosol layers fell due to gravity (Schoeberl et al., 2022). This amount of water vapor exceeded previous values of water vapor injected into the stratosphere, which is estimated at 10% of the actual quantity of water vapor amount in this layer (Millán et al., 2022). The large amount of water vapor converted ${SO}_2$ into sulphates and hydroxide (Sellitto et al., 2022), which reduced the proportion of ${SO}_2$ in the stratosphere (Zhu et al., 2022).

The amount of ${SO}_2$ which is injected in the stratosphere after the explosion was approximately 0.4Tg (terrogram) as measured by the TROPOspheric Monitoring Instrument (TROPOMI) (Zuo et al., 2022). (Millán et al., 2022) showed that this large amount of water vapor will remain for several years in this layer, which will directly affect the solar radiation budget and play a significant role in influencing the climate. However, this relatively small percentage of ${SO}_2$ also plays an important role on the climate global circulation (Zhu et al., 2022). It may have a greater impact if this volcano erupts again and injects more ${SO}_2$ into the stratosphere (Zuo et al., 2022). To complete these monitoring studies, new state-of-the-art Chemical Transport Models (CTMs) can estimate the volcanic ash's spatial distribution based on meteorological data inherited from global models at regular time steps covering the release time. We conducted simulations in this study using the WRF-Chem model (Grell et al., 2005; Fast et al., 2006) to investigate the causality of surface nudging effects on the meteorology and volcanic ash clouds results.

The volcano module within the WRF-Chem model (Stuefer et al., 2013) has shown good agreement with the available measurements like ground-based lidar as well as satellite data when conducted for typical volcanic eruptions (e.g. Mount Redoubt, 1989 and Eyjafjallajökull, 2010). Meteorological data assimilation is a process that uses observations of various weather variables to improve the accuracy of numerical weather prediction models.

Meteorological data assimilation is the process of integrating observational data with the beginning conditions of a model and consistently updating the model's state to enhance its depiction of the atmosphere. Different assimilation strategies, such as variational methods and ensemble-based methods, are used to address uncertainty in both the observations and the model. These strategies manipulate the model's variables, such as temperature, pressure, and humidity, in order to align them with the available data and generate meteorological forecasts that are more precise and dependable.

Data assimilation greatly enhances the precision of meteorological forecasts, leading to a wide range of practical uses. Precise weather forecasts in agriculture facilitate the optimization of planting, irrigation, and pest management techniques, resulting in enhanced crop productivity and minimized ecological consequences. Accurate predictions of wind, temperature, and precipitation are useful for planning routes, managing air traffic, and ensuring the safe functioning of road and rail networks in transportation. Moreover, in the realm of public safety, the provision of prompt and precise meteorological data is of utmost importance for the implementation of early alert systems, formulation of strategies for disaster response, and safeguarding of human lives and assets.

One technique for data assimilation is four-dimensional data assimilation (FDDA; Stauffer et al., 1990;1994), which continuously combines observational data and model results over a specific time step. FDDA is numerically inexpensive and can improve the accuracy of model outputs by adding artificial tendency terms to the model's prognostic equations and using suitable nudging coefficients to push the model towards meteorological observations.

The Flux-Adjusting Surface Data Assimilation System (FASDAS) is a relaxation mechanism that can be employed in the Weather Research and Forecasting (WRF) model to compel simulations to closely resemble a sequence of studies at every grid point. FASDAS utilizes the advanced data assimilation of using both direct and indirect methods to enhance the precision of boundary layer processes in numerical weather forecasts. FASDAS therefore modifies the heat, moisture, and momentum exchanges at the surface in the model using up-to-date observational data, thereby decreasing inaccuracies and enhancing the depiction of atmospheric phenomena.

FASDAS has proven to be quite advantageous in enhancing air quality modeling scenarios. FASDAS utilizes surface observations of contaminants such as ozone, carbon monoxide, and particulate matter to enhance the precision of air quality simulations. This capacity is essential for comprehending the influence of different emission sources, evaluating the efficiency of air pollution mitigation techniques, and delivering dependable air quality predictions.

The efficacy of the FASDAS relaxation method has been rigorously evaluated across a wide range of applications. Experts and scholars have assessed its efficacy in several meteorological and environmental situations, encompassing the prediction of severe weather occurrences, investigation of local climate trends, and evaluation of the dispersal of contaminants in urban regions(Alapaty et al., 2008; Gilliam et al., 2021; He et al., 2017; Osuri et al., 2020; Snoun et al., 2019b, 2021; Tong et al., 2020).The aforementioned studies have provided empirical evidence regarding the efficacy of FASDAS in augmenting model predictions and deepening comprehension of intricate atmospheric and environmental phenomena.

In its entirety, FASDAS signifies a substantial progression within the domain of atmospheric and environmental sciences. Through the integration of observational data and numerical models, FASDAS provides a robust instrument for enhancing atmospheric research, air quality modeling, and weather forecasting. Due to its well-established functionalities, it is an invaluable resource for meteorologists, air quality scientists, and researchers aiming to augment our comprehension of the Earth's atmosphere and its interactions with the surrounding environment.

Thus, in this study, the efficiency of meteorological data assimilation in improving the accuracy of atmospheric chemistry modeling was evaluated in the context of a volcanic eruption. The conventional datasets used to drive regional Numerical Weather Predictions (NWP) and CTMs are not always able to accurately estimate meteorological parameters for use in dispersion models. Thus, introducing more accurate meteorology through data assimilationcan improve the interactions between atmospheric processes, leading to a more realistic representation of model performance. This is the first mesoscale modeling study to examine this volcanic eruption.

The manuscript is organized as follows: a description of the submarine volcanic eruption on January 15, 2022 is provided in Section 2, the data assimilation methodology and model setup are discussed in Section 3. The results and discussions in Section 3 and Section 4 are divided into three subsections, which analyze the meteorological fields, volcanic ash and ${SO}_2$ concentrations based on WRF-Chem forecasts and observational data, and describe the evolution of the emitted plume using TROPOMI Sentinel5-P observations. Lastly, the discussions and conclusionsare summarized in the final section.

2. Materials and Methods

2.1. The volcanic activity

The Hunga Tonga-Hunga Ha'apai volcano (20.536°S; 175. 385°W) is located in the south Pacific Ocean and is 65 Km distant from Tongatapu, the main island of Tonga republic (Red triangle on Figure 3). Since 2014, it has remained relatively inactive till December 20th, 2021, when a multitude of volcanic eruptions has been recorded for the volcano. Interestingly, the large blast that rocked the volcano on January 15th and seen by the NOAA's Geostationary Operational Environmental Satellite 17 (GOES-17; Figure 1 and Supplementary Video S1), show the biggest explosion captured on camera. The eruption belongs to the category of submarine volcanic eruptions, which is a volcanic eruption that occurs under the surface of the ocean. These types of eruptions are often difficult to detect and can be dangerous due to the potential for tsunamis and other hazards. Submarine volcanic eruptions can also have significant impacts on marine life and the surrounding environment. Thus, it is important for scientists to monitor and study these types of eruptions to better understand their impacts and to be prepared for future events. It's evident that the submarine volcano violently exploded beginning from 04:30 (UTC) (Figure 1, left panel). The mushroom-shaped ash produced at 05:00 (UTC) sent pressure waves across the atmosphere (Figure 1, top middle panel), tsunami waves, and sonic booms (Lin et al., 2022;Yuen et al., 2022). After the plume expended radially at the top of the plume and reached its maximum extent (umbrella cloud of roughly 500km), the volcanic plume expanded upward and outward over the south Pacific, and began traveling westerly following the wind direction near 05:30 (UTC) (Figure 1, right panel). The plume topmost, as captured by satellite images, was at least 600 km in diameter by 19:03. The explosion was so huge that it obliterated the inhabited island of Hunga Tonga-Hunga Ha'apai (GVP, 2022).

2.2. Model and data assimilation descriptions

2.2.1. Flux adjusting surface data assimilation system (FASDAS)

(Alapaty et al., 2008) implemented the FASDAS surface and tropospheric nudging for mesoscale models to minimize the biases in meteorological predictions. The method continuously adjusts water vapor mixing ratio, surface air temperature, soil skin temperature and soil moisture in 3D space. FASDAS nudging is composed of a direct nudging of the surface atmospheric layer, and an indirect nudging via the soil moisture and temperature. For this purpose, in order to maintain thermodynamic stability within the model, the FASDAS approach adjusts atmospheric and land surface variables. In the FASDAS’s direct nudging approach, an error correction in the overall dynamical and physical processes that affects the surface moisture and temperature is performed. The indirect nudging, on the other hand, establishes adjustment of surface sensible heat fluxes to reduce the errors in soil moisture and soil temperature. Note that the FASDAS is only compatible with the Noah Land Surface Model (LSM); (Chen and Dudhia, 2001) and the Younsei University (Hong et al., 2006) Planetary Boundary Layer (PBL) schemes to run properly. Based on these facts, the equation for an adjusted FADSAS ground-skin temperature is given by:

$${\mathit{T}}_{\mathit{g}}^{\mathit{F}}={\mathit{T}}_{\mathit{g}}+\frac{\left({\mathit{H}}_{\mathit{S}}^{\mathit{F}}-{\mathit{\psi }}_{\mathit{q}}{\mathit{H}}_{\mathit{l}}^{\mathit{F}}\right)}{{\mathit{C}}_{\mathit{g}}}∆\mathit{t}$$

(1)

In Eq. (1), ${\mathit{H}}_{\mathit{l}}^{\mathit{F}}{\mathit{H}}_{\mathit{S}}^{\mathit{F}}$ are respectively the adjusted surface latent and sensible heat fluxes (W.${\mathit{m}}^{-2}$). The ${\mathit{T}}_{\mathit{g}}$and ${\mathit{T}}_{\mathit{g}}^{\mathit{F}}$ terms denote respectively the first guess and the assimilated soil temperature, ${\mathit{C}}_{\mathit{g}}$ represents the thermal capacity slice of the upmost soil per unit area (J.${\mathit{k}\mathit{g}}^{-1}$.${\mathit{K}}^{-1}$), $∆\mathit{t}$ is the time step and ${\mathit{\psi }}_{\mathit{q}}$ represents a normalized weighting factor for soil moisture adjustment.

2.2.2. Model scenarios configuration

The inline coupled WRF-Chem mesoscale model v4.4 was used to provide meteorological and volcanic ash predictions for regional applications on a countrywide scale or smaller. The compilation uses the self-installation scripts provided by (Hatheway et al., 2023).For this study, the model covered a single computational domain with 88 km horizontal resolution covering the western Pacific, Tonga, New Zealand, and Australia (Figure 2). The grid size for the domain was 95×50 points (E-W by N-S). The ARW core v4.4 configuration used a hybrid vertical coordinate (Park et al., 2013) which quickly removes the influence of terrain on higher coordinate surfaces while remaining terrain following (TF) near the surface (Snoun et al., 2019a). The projection method used was the Lambert Conic Conformal (LCC). The meteorological forecasts were initialized with hourly ERA5 lateral boundary data at 0.25°× 0.25° horizontal resolution and hourly time step. Hourly WRF-Chem ground level meteorological outputs were collected at 4 surface stations (shown as yellow dots in Figure 2, Table 1) for comparison. To investigate the impact of data assimilation on predicted meteorology and volcanic ash, three WRF-Chem runs were conducted: a base case scenario (OMET) which only considered meteorology and the volcanic eruption, a second scenario (NODA) which included the volcanic eruption and Goddard Chemistry Aerosol Radiation and Transport (GOCART) background data but no data assimilation, and a third scenario (FASDAS) which included nudging, gas andaerosol chemistry using natural, biogenic, and anthropogenic emissions in addition to the volcanic eruption report.

Major physical parameterizations settings involve the WRF double moment Morrison microphysics scheme (Morrison et al., 2009), the surface physics option from the NOAH Land Surface Model (LSM; Chen and Dudhia, 2001), the revised similarity theory of Monin-Obukhov for surface layer (Jiménez et al., 2012), the Rapid Radiative Transfer Model (RTTMG) parameterization for shortwave and longwave radiations (Iacono et al., 2008),the Grell-3D cumulus scheme (Grell and Dévényi, 2002) and the Yonsei University (YSU) planetary boundary layer (Hong et al., 2006). To model volcanic ash emissions from the Hunga Tonga volcano, we used the pollutant emissions numerical tool prep_chem. This software provides information on the location and timing of the eruption (January 15th, 16:15 UTC) and ash height. Based on research by NASA scientists and (Proud et al., 2022), we used a 58 km ash plume height. The forecasts also employed a simple dust treatment that considers ash fall and concentrations distributed over 10 particle size bins (chem_opt=402) to better simulate the transport and fall of ash of different sizes, and the S2 classification (Stuefer et al., 2013; Table 2).

We conducted a 15-day forecast from January 10th to 25th, 2022 on an hourly basis for both FASDAS and non-data assimilated runs, with a 5-day model spinup period to improve consistency in model physics and dynamics. The model includes 45 pressure levels, with the top level at 5000 Pa. The FASDAS system (Alapaty et al., 2008) is used to simulate boundary layer nudging. We choose relaxation terms of 8.3×${10}^{-4}{\mathit{s}}^{-1}$for the surface, water vapour and temperature nudging. The radius of influence is set to be 2×Δx (i.e.176 km), following the recommendations of the FASDAS user guide. In addition, according to (Bullock et al., 2014) sensitivity tests, the free atmosphere relaxation coefficients are selected as 5×${10}^{-6}{\mathit{s}}^{-1}$for moisture, and 5×${10}^{-5}{\mathit{s}}^{-1}$ for temperature as well as u and v wind components. To improve the performance of the model, FASDAS nudging input analysis data was adjusted to a higher temporal frequency of hourly intervals and fed to OBSGRID. This approach, previously discussed in (Li et al., 2016), has proven to be more effective than the default 3 hourly frequency of input analyses. We also performed objective analysis of surface pressure (Reen et al., 2015) to prevent the WRF model from rejecting the surface analysis as being too distant from the actual surface, and to improve the estimation of the surface pressure. In addition, we used the approach of (Reen et al., 2016) to decrease the magnitude of drying on each Cressman scan analysis, reducing overdrying when the diagnosed drying need at a specific location is spread to a second point where the first guess is drier than the first guess at the observation's location. The contours were partly generated using the PostWRFopen-source toolkit (Nikfal, 2023) for easy visualization.

Additionally, we computed six model performance statistical metrics for model evaluations (with and without assimilation) including the FAC2 (fraction of model in a factor of two of actual observations), the RMSE (root mean squared error), the MB (mean bias), the normalized mean bias (NMB), the index of agreement (IOA) and Pearson correlation (r). The metric indexes are expressed as:

$$\mathrm{M}\mathrm{B}=\frac{1}{\mathrm{n}}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{M}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}$$

(2)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}={\left[\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{M}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right)²}{\mathrm{n}}\right]}^{\frac{1}{2}}$$

(3)

$$\mathrm{r}=\frac{1}{\mathrm{n}-1}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left(\frac{{\mathrm{M}}_{\mathrm{i}}-\overline{\mathrm{M}}}{{\mathsf{\sigma }}_{\mathrm{M}}}\right)\left(\frac{{\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}}{{\mathsf{\sigma }}_{\mathrm{O}}}\right)$$

(4)

$$\mathrm{FAC}2=\mathrm{data}\mathrm{portion}\mathrm{verifying}0.5\le \frac{{\mathrm{M}}_{\mathrm{i}}}{{\mathrm{O}}_{\mathrm{i}}}\le 2$$

(5)

$$\mathrm{I}\mathrm{O}\mathrm{A}=\left\{\begin{array}{c}1-\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{M}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right|}{\mathrm{c}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right|},\mathrm{if}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{M}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right|\le c{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right|\\ \frac{\mathrm{c}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right|}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{M}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right|}-1,\mathrm{if}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{M}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right|>c{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right|\end{array}\right.$$

(6)

$$\mathrm{N}\mathrm{M}\mathrm{B}=\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{M}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{O}}_{\mathrm{i}}}$$

(7)

The above simulated and observed values are denoted as ${\mathit{M}}_{\mathit{i}}$ and ${\mathit{O}}_{\mathit{i}}$, respectively.

3. Results

3.1. Impact of Data Assimilation on Meteorological Fields

There is a potential connection between atmospheric depression and underwater volcanic eruptions. When an underwater volcanic eruption occurs, it can release large amounts of gas and ash into the atmosphere. This can lead to the formation of a low-pressure system, or atmospheric depression, as the gases and ash can block out sunlight and create cooling effects in the surrounding area. Additionally, the release of gases and ash can also disrupt wind patterns and cause changes in atmospheric pressure. This can lead to the development of storms and other weather phenomena, which can further contribute to atmospheric depression. Thus, it seems that there may be an indirect influence of this underwater volcanic eruption on atmospheric pressure systems. Analysis of pressure systems during 4 days including the eruption event shows that it reinforced a disturbance of the pressure system, where a band of low pressure is apparent to the east in the vicinity of the volcano in all the plots of Figure 3. In fact, sea level pressure of WRF-Chem and ECMWF's Fifth generation Reanalysis (ERA5) data analyzed from 14-17 January 2022 (Figure 3a–h) show that mesoscale forecasts of wind flow and pressure patterns suggest the presence of prevailing tropical westerly wave motion, which is combined by the movement of two air masses: (i) A tropical warmer air mass moving from higher latitudes (1020-1036 hPa) and (ii) cold polar air masses coming from lower latitudes (990-1020 hPa). Under the influence of this situation, the flow field shows anticyclonic vorticity localized in New Zealand and the south-western part of Australia (Figure 3a–d), which fosters southerly circulations in the domain. Associated with these vorticities, movement and formation of cyclonic (i.e. low pressure) systems are noted in the northern part of Australia and New Zealand, in agreement with the data reported in the ERA5 reanalysis (Figure 3e–h). Under the influence of localized low pressures, the 990 hPa flow field on 14-17 January produced north-westerly winds nearby the Hunga Tonga volcano, and a noticed alleviation of the low-pressure systems the days following the eruption.

The ERA5 data (Figure 3e–h) show similar pressure levels with those produced with WRF-Chem. This demonstrates that the flow field with the regional model does not practice a topographic influence, and consequently follows the ERA5 circulation. Thus, the formation of the storm (the well dug depression of 14-17 January) makes it relevant the link between the large amounts of ash and particles into the atmosphere, which can affect weather patterns and cause formation of low-pressure systems during the large eruption of the Hunga Tonga volcano. However, the exact impact of an underwater volcanic eruption on atmospheric pressure systems needs further investigation, as it depends on a variety of factors, including the size and intensity of the eruption, the location of the volcano, and the atmospheric conditions at the time of the eruption.

Figure 3. Simulated WRF-Chem (left panels) sea level pressure (hPa) and data reported from ERA5-reanalysis (right panels) during 14-17 January 2022.

Figure 3. Simulated WRF-Chem (left panels) sea level pressure (hPa) and data reported from ERA5-reanalysis (right panels) during 14-17 January 2022.

Meteorology drastically influences the formation and transport of atmospheric aerosols. We selected three meteorological parameters for model evaluation: temperature, wind speed and direction, and four random surface stations (Fuamutu: 21.14S; 175.09W, La Tontouta: 22.01S; 166.13E, Auckland: 37.01S; 174.48E and Brisbane: 27.23S; 153.06E). Figure 4 depicts the time series of measured and model-retrieved hourly 2 m temperature (°C). The statistical comparisons summary is shown in Table 3. The model reproduces quite well temporal variations of temperature, with the FASDAS run showing more accurate results through assimilation of observational data. For instance, compared to OMET and NODA cases, we achieved 56.25% and 55% reduction in MB, and 28% and 32% improvement in IOA, respectively. The FASDAS also accurately reproduces temperature with 95% correlation, resulting in enhancements of 10.46% and 18.75% over OMET and NODA, respectively.

Figure 5 examines the comparisons of wind speeds and directions at 10m for the three meteorological datasets. The OMET, NODA and FASDAS winds are positively biased compared with the wind roses reported from the measurements averaged over the surface stations. The figure also clearly shows the mean positive bias in wind directions (except for FASDAS) following the northern direction and not exceeding 15%, and a higher proportion of positively biased wind speeds between 0-5 ${m.s}^{-1}$. The FASDAS biases for wind speeds are smaller (NMB= 29%) as compared to OMET and NODA (NMB= 33 and 32%, respectively). Thus, the WRF-Chem predictions systematically overestimated the wind speed by 31% on average. This bias probably results from the domain's coarse resolution (Yahya et al., 2015; Zhang et al., 2015). The evolution of atmospheric volcanic ash particles and aerosols is sorely affected by meteorological fields such as boundary layer height, incoming solar radiation and air temperature. Figure 6, Figure 7 and Figure 8 illustrate the mean impact of data assimilation on Planetary boundary layer height, 2 m temperature and downward shortwave flux at the ground over Oceania and Australia in January 2022. Figure 6 shows that inclusion of the chemistry mechanism in the simulations (case NODA) was found to only have minor influence on the predictions of mean PBL Height. The FASDAS effect was more pronounced than NODA and OMET cases, and produced 300-800 m lower PBL Height over central Australia and the western part of Fiji Islands (Figure 6). Higher mixed layer depths with FASDAS were predicted for some areas around the erupted volcano, over Australia and New Zealand. Over the Australian extreme easterly and westerly land portions, the simulated increases in mean PBL Height for the FASDAS case are associated with 100-200m higher atmospheric boundary layer height (Figure 6, middle and bottom panels), which could be attributed to the prevalence of relatively cooler air masses above the pacific sea. The changes in mixed layer depths also affect the 2m temperature rates. Average temperatures simulated in the western part of the model domain were till 6°C lower for the FASDAS case than for the only meteorology (OMET) and baseline (NODA) cases, respectively (Figure 7). Higher temperatures (up to 5 °C) were predicted for most regions in the eastern part of Australia and west of 120°E for the FASDAS case. The change in the temperature distribution is partly reflected in the abundance of the surface temperature measurements assimilated during the FASDAS process, plus thespatial distribution of the mean simulated atmospheric boundary layer height (Figure 6) that displays a sharp decrease over central Australia for the FASDAS case. In general, as the temperature of the surface increases, the PBL will become more unstable and will tend to rise. This is because warmer temperatures increase the amount of heat being transferred from the surface to the atmosphere, which in turn increases the amount of turbulent mixing. Conversely, when the temperature of the surface decreases, the PBL will become more stable and will tend to descend. Thus, temperature is an important factor that can affect the height of the PBL, and that the large differences in the results may indicate the effectiveness of the FASDAS effect in retrieving more accurate temperature (and consequently PBL height) than the without assimilation cases.

Similarly, the inclusion of FASDAS assimilation (case FASDAS) was found to have significant influence on the 10-day averages of solar radiation. Figure 8 depicts the 15-25 January mean incoming solar radiation, including nighttime periods. For regions with high volcanic ash emissions (around the volcano as well as northern New Zealand), simulations for cases without assimilation predicted less solar radiationrates than the FASDAS case (50-150${W.m}^{-2}$), whilst the opposite case is observed at the locations where temperature and PBL height positive anomalies have been detected, which prove that the temperature and PBLH distributions are drastically reflected in the pattern changes of the solar radiation. In fact, when the amount of shortwave radiation (or sunlight) increases, the temperature of an object or area also increases. This is because the shortwave radiation is absorbed by the object or area and converts into heat, causing the temperature to rise. Conversely, when the amount of shortwave radiation decreases, the temperature also decreases (the greenhouse effect), in which shortwave radiation is absorbed by the Earth's atmosphere and surface, leading to an increase in temperature, which is rigorously reflected in the contours of Figs. 7-8.

1.2. FASDAS effect on volcanic ash and ${SO}_2$ concentrations

Volcanic ash concentration can significantly influence meteorological data assimilation because it can affect the accuracy of observations and measurements taken by various sources. For example, ash particles in the atmosphere can interfere with satellite sensors and affect the accuracy of satellite-based observations. Ash particles can also affect the accuracy of measurements taken by weather stations and aircraft. In addition, volcanic ash can have a significant impact on the weather itself, as it can block out sunlight and reduce solar radiation reaching the Earth's surface. This can result in lower temperatures and increased precipitation in affected areas. To quantify the effectiveness of the FASDAS effects, Figure 9 presents the comparison between the total ash columns transport of ash particles. Output from NODA and FASDAS scenarios in WRF-Chem, specifically the finest three ash bins in the ${PM}_{10}$size range. It's important to note that the calculation of ${PM}_{10}$ when incorporating volcanic ash scheme involves aggregating the quantities of the smallest ash particles from the distribution referenced in Table 2, specifically those particles with a diameter smaller than 10 micrometers. The formula for calculating ${PM}_{10}$ is based on the total of these fine ash particles:

$${PM}_{10}(\mu g.m^{-2})=\left[\frac{1}{2}\times vas{h}_8+vas{h}_9+\frac{1}{2}\times vas{h}_{10}\right]\times {\rho }_{air}$$

(8)

The Volcanic ash bins are measured in micrograms per kilogram (unit $\mu g.{kg}^{-1}$), and the density of air, denoted as ${\rho }_{air}$, in kilograms per cubic meter ($kg.m^{-3}$). The outcome, the derived surface concentration of ${PM}_{10}$, is then presented in micrograms per square meter ($\mu g.m^{-2}$). This measurement effectively represents the concentration of ash across a column of air.

As depicted in Figure 9, the transport fields of ash cloud particles reconstructed by NODA and FASDAS runs are consistent for some parts of the domain, with an apparent more predicted ash with nudging. Nearby the volcanic eruption location, the volcanic ash predictions are mostly associated with rate differences. The NODA experiment generally underestimated the ash retrieved by FASDAS, though it shows similar spatial orientation. Thisunderestimation tendency is caused by the lower concentration plume predicted without data assimilation, which attracts the predicted ash particles to accumulate in this area and favors anegative bias. Another factor that could be associated with the differences is the simulation of lower wind speed and PBLH, which limits the movement of air masses and promotes the formation of pollution peaks.

In summary, the difference between the results of NODA and FASDAS scenarios is rather at the level of the plume itself. For our situation, the plume was more elongated due to atmospheric instability limiting dispersion, which gave a faster propagation due to a higher wind speed, as ash will disperse more widely and for longer periods of time in unstable atmospheric conditions, such as those associated with strong winds and low-pressure systems. Therefore, it is important to consider the concentration of volcanic ash when performing meteorological data assimilation in order to improve the accuracy of weather forecasts and to better understand the impact of volcanic ash on the weather.

We assessed the precision of ${SO}_2$ surface levels both before and after data assimilation through the use of local ${SO}_2$surface observations. Monitoring stations, mainly from Australia, were evaluated during the simulated period as they provide the most dependable data. A comparison of the spatial distribution of ${SO}_2$ surface concentrations against observations were presented in Figure 10, it can be observed that NODA has a less concentrated ${SO}_2$ cloud than FASDAS (as seen in Figure 10 top panel) but both models predominantly underestimated the observed locations.

The results of the simulation showed that neither NODA nor FASDAS displayed significant ash concentrations at the evaluated locations. It was found that after incorporating ground site measurements, the concentration of ${SO}_2$ was significantly increased and the ash cloud appeared to be spreading further to reach the Fiji islands, indicating that the model is utilizing more information about the atmospheric conditions that influence the ash cloud, specifically an increase in ash production (as seen in Figure 10-bottom panel). However, it was also noted that the ${SO}_2$ concentrations were still underestimated in areas with limited site observations, such as at the Gap station (18.6 ppb observed against 5 ppb predicted). This underestimation is likely due to unrealistic model anthropogenic emissions, PBL height estimations and more importantly errors due to ash cloud retrieval within the chemical mechanism employed in the model (further discussed in the conclusion).

In fact, as the boundary layer height estimation decreased in FASDAS (following Figure 6), the concentration of ${SO}_2$ near the surface may increase because the ${SO}_2$ has less space to disperse. Indeed, on a FASDAS cooler (in agreement with Figure 7) and humid day, the boundary layer may be lower due to the reduced convective mixing. It actsas a mixing layer, bringing together different air masses and promoting the exchange of gases.

1.3. Analysis of Sulfur dioxides transport using TROPOMI/S5P

To understand the state of the plume of volcanic ash from the days following the explosion, a spatial analysis of TROPOMI/S5P standard level 2 (L2) ${SO}_2$ column patterns ($DU$) beingpresented. This is important because in addition to the volcanic ash emitted, a large quantity of ${SO}_2$ is released into the atmosphere at the time of the eruption. The input dataset for this analysis is mainly the high-resolution record of the instrument data acquired over the study period of 14-21 January 2022, at a daily revisiting interval. The investigated areas include the whole Tonga archipelago and Australia for the sulfur dioxide plume. The panel plots of Figure 11 summarize the satellite's daily measurements available for cross-comparison.

Because of the less significant sub-aerial eruption that occurred on January 14, 2022 (Figure 11a), there is evidence of ${SO}_2$ volcanic plumes in the Northern direction of the volcano. In fact, before the major eruption of January 15th, a paroxysm with a slightly lower amplitude occurred the previous day (January 14, 2022). The submarine larger eruption that occurred the following day (Figure 11b) produced expanding eruption plume that started to cover nearby Tongan islands with ${SO}_2$. The apparent rates were supposed to be greater than that plotted in the figure, but the amounts do not appear to be as extensive as expected. This is due to the moment of the sensing time (00:05 am, UTC) of the TROPOMI/S5P on January, 15th that preceded the detected massive blast. Indeed, ${SO}_2$ cloud ejections became more apparent on January, 16th (Figure 11c) hitting the islands of Fiji, Vanuatu, and New Caledonia, and the thresholds reported in the plume center reached 15 DU.

Monitored column concentration of sulfur dioxides in the following days shows that the plume was moving with the westerly wind direction (Figure 11d) and penetrating the Queensland region across the far northeast of Australia, before reaching Australia's Northern Territory, Queensland, and Western Australia on January 18th (Figure 11e). In particular, on January 17th, the intense ${SO}_2$ plume center which lies in the pacific sea near the east coast of Australia continued to spread westward, resulting in crossing the coast of Australia until it is transported towards the above-mentioned locations in the country. During the period 19-20 January, the maps report a lightning in the intensity of the plume and continuation of the propagation in the westerly direction (Figs. 11f-g). The plume has still covered north and northwestern parts of Australia, whilst it entirely shifted in destination to the east and central Java parts of the Indonesian country (Figure 11h). The recorded plume column density values are fluctuating between 1-4 DU. Finally, because there is no record of extra eruptions in the days following January 19th, the atmospheric circulation has played the role of dispersing the ${SO}_2$ plume while attenuating its intensity over space and time.

4. Discussion

Volcanic ash emissions can remarkably impact the environment and the air we breathe at different atmospheric scales, depending on the event's extent. The Hunga Tonga-Hunga Ha'apai volcanic eruption was so violent that caused a tsunami that spread across the Pacific Ocean and whose shock wave circled the earth several times. Through the application of the WRF-Chem model, we have depicted a systematic study of the Hunga Tonga-Hunga Ha'apai eruption recorded in the middle of January 2022, along with its role in volcanic ash and pollutants transport in the pacific area. In addition, since meteorological parameters (and consequently pollutant concentrations) can change in distribution because of the presence of aerosols in the atmosphere, we investigated in this study an implemented FASDAS parameterization suitable for WRF-Chem v4.4 coarser domains. An ultimate objective from the experiments was to effectively comprehend the direct effects that can exert the usage of FASDAS nudging on near-surface wind patterns, temperature and solar radiation.

The inclusion of FASDAS in the chemistry runs results in a pronounced alleviation of PBL heights (up to 800 m), 2m temperature (up to 6 °C) and solar radiation (up to 250 W.$m^{-2}$) profiles. These reductions were associated with a higher mean ${SO}_2$ and volcanic ash concentrations. Evaluation of model results with surface weather station measurements gave slightly more biased results in wind speed with NODA (overall NMB=32%) than FASDAS (overall NMB=29%). WRF-Chem produced negative biases for temperature (NMB=4% with NODA against 2% with nudging). The biases reflected the importance of the impaction of nudging effects. Furthermore, the comparison of the volcanic ash cloud diffusion with and without assimilation revealed differences near the eruption site, mainly because of the predominant wind speed overestimation rates by the two model runs.

The vash biases (not shown) suggest the role of surface nudging effect in the spreading out of Hunga Tonga-Hunga Ha'apai volcanic ash until the New Caledonia islands and Australia. The study also encompasses a spatial sampling analysis step, which is aimed to follow the plume's movement in the days following the eruption. Our study, with the help of TROPOMI/S5P ${SO}_2$ column data, suggested that the volcanic ash plume drifted from the main eruption site and the circulation underwent the Western direction to reach Australia, which comply with the wind direction simulated by the model. The Hunga Tonga-Hunga Ha'apai volcanic eruption made a notable impact on the Earth's atmospheric composition, particularly concerning the water vapor mixing ratio. According to the research published by (Xu et al., 2022), the eruption caused a considerable injection of water vapor into the stratosphere, which is unusual for volcanic events, as they typically emit large amounts of sulfur dioxide (SO2) but only minimal water vapor. In contrast, the Hunga Tonga eruption was distinctive due to its minimal SO2 and significant water vapor emissions. This unique characteristic led to an increased concentration of water vapor in the stratosphere, which has implications for atmospheric chemistry and potentially climate patterns.

The findings from (Schoeberl et al., 2022) further corroborate these observations, emphasizing the rarity and significance of such a high water vapor content being introduced to the stratosphere from a volcanic eruption. This event, therefore, stands as a notable exception in the study of volcanic eruptions and their impacts on the Earth's atmosphere. Figure 12 represents a statistical analysis of the water vapor content differences between two distinct dispersion simulations of the underwater eruption. The mean value (μ) of the water vapor content difference is approximately 3.52 x ${10}^{-5}$kg ${Kg}^{-1}$, which is quite small, suggesting that the overall impact of data assimilation on the modeled water vapor content is minimal but detectable. The standard deviation (σ) of approximately 0.00093 kg ${Kg}^{-1}$ indicates a relatively tight clustering of differences around the mean, which reflects a consistent modification of the water vapor field by the data assimilation process. Additionally, the histogram’s bell-shaped curve and its symmetry about the mean suggest a normal distribution of differences. This pattern suggests that the data assimilation process used in the FASDAS simulation has introduced systematic adjustments to the water vapor content, which could be due to the assimilation of observational data improving the representation of the volcanic plume and its dispersion in the atmosphere.

5. Conclusions

In light of recent studies, it is important to acknowledge the limitations of our current study regarding the modeling of volcanic ash dispersion from the Hunga Tonga-Hunga Ha'apai eruption. Notably, several recent publications, including (Sellitto et al., 2022, Legras et al., 2022, and Zhu et al., 2022), have presented evidence suggesting rapid sedimentation of ash and the dominance of sulphate aerosols, which indicates that the volcanic ash may not have dispersed significantly but remained localized. This challenges the foundational assumption of our study that there was substantial volcanic ash dispersion.Additionally, it appears that the WRF-Chem model used in our study may not adequately account for the rapid sedimentation of ash, potentially influenced by ice particle formation due to a significant water vapor injection during the event. This factor, along with possible inaccuracies in ash size distribution within the model, could have led to an underestimation of ash dispersion in our simulations, considering that the ${SO}_2$ and ash plumes are not co-located as previously assumed.This new understanding underscores the need for further refinement in modeling techniques and the incorporation of more accurate ash sedimentation and distribution parameters. Future studies should focus on improving the accuracy of volcanic ash dispersion models by considering the rapid sedimentation and distinct characteristics of ash and sulphate aerosols, as well as the unique conditions of each volcanic event.