High–Resolution Modeling of Airﬂows and Particle Deposition over Complex Terrain at Sakurajima Volcano

: The realistic representation of atmospheric pollutant dispersal over areas of complex topography presents a challenging application for meteorological models. Here, we present results from high–resolution atmospheric modeling in order to gain insight into local processes that can affect ash transport and deposition. The nested Weather Research and Forecasting (WRF) model with the ﬁnest resolution of 50 m was used to simulate atmospheric ﬂow over the complex topography of Sakurajima volcano, Japan, for two volcanic eruption cases. The simulated airﬂow results were shown to compare well against surface observations. As a preliminary application, idealized trajectory modeling for the two cases revealed that accounting for local circulations can signiﬁcantly impact volcanic ash deposition leading to a total fall velocity up to 2–3 times the particle’s terminal velocity depending on the size. Such a modiﬁcation of the estimated particle settling velocity over areas with complex topography can be used to parametrize the impact of orographic effects in dispersal models, in order to improve ﬁdelity.


Introduction
The transport and deposition of atmospheric pollutants over regions of complex topography is dictated by local-scale circulations enforced by the terrain, turbulent mixing within the planetary boundary layer (PBL), and the physical properties of the pollutants [1][2][3][4]. The successful reproduction of observed patterns hinges on the correct estimation of the involved particles' fall velocities [5]. This can be especially challenging in the case of volcanic ash, as it is introduced in a range of particle sizes (from 2 mm down to a few microns) with disparate particle characteristics (i.e., density, sphericity), and a large range of particle terminal velocities O(10 -2 -10 1 m s -1 ).
The simulation of volcanic ash transport and deposition (VATD) is especially challenging in the proximal region to the vent (∼5 km) owing to a combination of unresolved volcanogenic effects [6][7][8] and the impact of the volcano's orography. An orographic flow is the result of the atmosphere impinging on a topographic maximum and critically depends on atmospheric stability and the kinetic energy of the incoming flow. Under stable conditions, internal gravity waves (or mountain waves) [9] and strong downslope winds [10] can be triggered on the lee-side. These involve vertical velocities comparable to the horizontal wind speed and can play a vital role in ash deposition by modifying the particle fall velocities. However, the overall deposition distances typically studied (∼100 km) [11,12] only allow for computational modeling using relatively coarse resolutions (horizontal grid spacing, ∆x, over 2 km [13][14][15][16][17]), leading to a negligible resolved vertical velocity range and particle fall velocities approximately equal to the particle's terminal velocity along the whole deposition path.
Recently, a number of studies employing mesoscale modeling have shown that at least partially resolving orographic effects can be important for VATD modeling [18,19], but no attempt was made to quantify the effect. Additionally, the use of horizontal grid spacings within the turbulence gray area [20], i.e., the resolution range in which the dominant length of the flow is comparable to resolution, typically between 100-1000 m [21], was a limiting factor, as it comes with the risk of introducing convection generated by numerical noise [22] or numerical artefacts such as unrealistic oscillations in the output [23,24].
Fully resolving orographic effects from isolated mountains necessitates a large-eddy simulation (LES) approach, which allows for the direct resolution of the largest energycarrying eddies. This approach has been successfully used to resolve the detailed forcing of complex topography [25,26] or urban terrain [27][28][29][30][31] and to study the transport of atmospheric pollutants [32,33]. LES modeling applied to real case studies commonly uses mesoscale models for initialization profiles [33,34], but most LES studies include some degree of idealization; limiting the atmospheric structure, assuming horizontal homogeneity, and using periodic lateral boundary conditions.
In recent years, efforts are made to bridge a gap between mesoscale and LES modeling. Using LES models as the basis, this involves using increasingly realistic conditions such as a temporally varying forcing [35][36][37]. An alternative approach is multiscale modeling, which has been used to bridge the synoptic and local scales [24,34,38,39]. In multiscale modeling, a coarse LES is carried out concurrent to mesoscale modeling, skipping gray zone resolutions. The Weather Research and Forecasting (WRF) model [40] has been used in this way in a number of cases, including the study of orographic processes over complex topography [24,[41][42][43][44] and pollutant transport [2]. Despite the relatively coarse resolution compared to traditional LES modeling, multiscale modeling provides a flexible approach to study real cases without the high degree of idealization.
Here, the WRF model is employed in a multiscale setting to study the airflow and particle movement in the planetary boundary layer (PBL) over the area surrounding the Aira caldera and Sakurajima volcano, located on the island of Kyushu, southern Japan ( Figure 1). In this area, multiple peaks form an elongated ridge along the south-north direction, with a half-width of approximately 2 km and a maximum height of 1117 m and with a number of ridges formed along the western side, creating a challenging topography for numerical modeling. Current activity from the volcano is characterized by small plume heights (typically 1-3 km above ground level (agl) [19]), that are strongly influenced by the PBL structure (for example see Figure 2a) and orographic effects ( Figure 2b). Thus, the aim of the study is two-fold: (i) to assess the WRF model's capability of providing realistic simulations over the volcano, and (ii) to explore the possible impact of the resolved flow on the fall velocity of the ash particles.

General Description of Case Study Days
The study focuses on the early morning hours of two days marked with volcanic eruptions: 6 June 2017 and 1 October 2017. The June eruption occurred at 0756 Japan Standard Time (JST; JST = UTC + 9 h), with a short initial phase (1-3 min; ash ejected over 3 km and concealed by cloud cover), followed by intermittent emissions a few hundred meters agl over the following hour. Ash from the initial phase was transported eastwards, while ash from the second phase was transported westwards, owing to the synoptic conditions ( Figure 3a). On the day of the eruption, the Japan Meteorological Agency (JMA) reported that Kyushu entered the rainy season, a weather pattern dominated by the Baiu front; a stationary warm front south-southeast of the Japanese islands caused by the interaction of near-surface cool polar masses and warm tropical masses from the southwest. This is evident in the vertical structure of the atmosphere (Figure 3b). Close to the surface, wind direction is northeasterly, with the wind turning geostrophic at approximately 900 hPa. The two air masses are mixed between 800-700 hPa, leading to a layer of significant directional shear, which is capped with a temperature inversion and a cloud layer at 700 hPa. Above the inversion, wind direction is westerly, and the moisture content is increased. The October eruption occurred at 0613 JST and lasted for ∼1 h, featuring a volcanic plume height of ∼1 km agl that was transported towards the west, with the effect of vertical wind shear apparent due to the significant fanning of the ash cloud ( Figure 2a). The synoptic conditions were very similar to the June case, with circulation dictated by a combination of a high-pressure system to the east and a low-pressure system to the west (Figure 3c). Close to the surface observed wind directions were northeasterly, with the wind turning geostrophic at approximately 800 hPa, accompanied by a temperature inversion. Unlike the June case, wind direction was southerly with height, turning westerly above 600 hPa owing to the influence of the subtropical jet stream [45].
The vertical structure for both days is captured well in the European Reanalysis 5 (ERA5) [46] dataset (shown in black in Figure 3b,d). Despite the overall match between the ERA5 and sounding data, significant disagreement can be seen close to the surface (with the ERA5 data failing to reproduce the northeasterly surface flow) and below the change to the westerly flow (∼ 800 hPa in Figure 3b and ∼ 600 hPa in Figure 3d). Other than wind direction, in the October case, the ERA5 data overestimate the temperature lapse rate close to the surface. For both days, the wind shear structure is consistent between 0400-0900JST. Up to 1 km above sea level (asl), the wind shear profile for both days is very similar, albeit with a ∼15 • difference ( Figure 4). Differences between the two profiles are evident above 1 km asl as the sharp turn in wind direction is evident in the June case ( Figure 4a) while in the October case, wind direction stays between 110-150 • (Figure 4b). On both days, a Doppler lidar was placed on the west (lee with respect to low-level flow) side of the volcano (see Figure 1b and Section 3.1). For both cases, the perturbation part of vertical velocity (based on a 30 min moving average) revealed oscillations captured between ∼100-300 m agl ( Figure 5). During the June case, these oscillations were mainly concentrated around a period (T) of 10 min (Figure 5a), while in the October case, the main peak in amplitude is associated T = 3 min (Figure 5c). Waveforms presumed to be associated with gravity waves can be clearly seen in the signal (Figure 5b,d) after a moving average of T/3 (i.e., 3 and 1 min, respectively) has been applied. At an average advection speed 3.4 and 7 m s -1 for the June and October cases, respectively, these values lead to a length scale of 2040 and 1260 m (≈O (1 km)), which are comparable to the volcano's geometrical characteristics.

Observational Equipment
Acquisition of in-situ data over mountainous terrain is always a challenge [47] and volcanic activity further complicates the issue [48]. Owing to the long-lasting activity of the Sakurajima volcano (since 1955) and the health concerns for the residents in the vicinity [49], an automated observational network is operated near the volcano. In order to minimize ambiguity in the atmospheric state rawinsonde data, surface station and averaged Doppler lidar data were used to characterize the stationary response of the flow, while instantaneous Doppler lidar data are used to assess instantaneous winds and the magnitude of turbulence.
Surface wind speed and direction are measured by two nation-wide systems (circle markers in Figure 1b): seven stations from the Ministry of Environment's (MoE) Atmospheric Environmental Regional Observation System (AEROS or Soramame in Japanese, http://soramame.taiki.go.jp/ (accessed date 25 February 2021); see [3] for details) and three stations from the Automated Meteorological Data Acquisition System (AMeDAS, https://www.jma.go.jp/jma/en/Activities/amedas/amedas.html (accessed date 25 February 2021)). All data are hourly, but AEROS data are instantaneous while AMeDAS data are 10-min averages of the values before the observation. Doppler lidars have been commonly used to measure wind speed and turbulence characteristics in a number of cases in mountainous terrain [50,51] and in harsh environments [52,53]. An optical-fiber Doppler lidar (LR-07F type III S; Mitsubishi Electric Corporation, Japan) is continuously operating on Sakurajima [54], located near the top of Mt. Haruta (31.59 • N/130.63 • E; star marker in Figure 1b). The wavelength is 1.5-1.6 µm, and the lidar is set pointing upwards, with a vertical resolution of 30 m between 30-600 m agl. Usable data usually lie within a thinner layer between 100-300 m agl (∼500-700 m asl). The laser pulse width is 200 ns and the pulse repetition rate is 4 kHz. This leads to an effective frequency of ∼0.1667 Hz (i.e., ∼1 measurement per 6 s). Here, a 30-sec moving average was applied to remove noise from the data and values were then resampled to one observation per minute (∼0.0167 Hz).

WRF Model Settings
Simulations were carried out using WRF version 4 [40]. As noted in Section 2, the dominant length scale for a stable PBL is unresolvable here. Instead, the horizontal grid spacings for the LES domains were chosen to represent gravity wave related instabilities, with dominant length scales of a few hundred meters [21]. Four one-way nested domains were used in total, progressively focusing over the volcano ( Figure 6). For the two outermost domains, D1 and D2 (Figure 6a), with both 166 by 166 horizontal grid points at ∆x = 3150 and 1050 m, respectively, a PBL parametrization scheme [55] was utilized, while for the two innermost domains, D3 and D4 (Figure 6a), with 421 by 421 horizontal grid points at ∆x = 150 and 721 by 394 horizontal grid points at ∆x = 50 m, respectively, an LES scheme [56] was used. No turbulence-seeding technique [34,39] was used, as turbulence seeded at the boundaries is expected to diminish with distance. The reason why one-way nesting is used is because the horizontal sizes of the inner domains are limited and, thus, any resolved features in the inner domains will be quickly wiped out by the incoming flows at the lateral boundaries or the effects of the feedback from an inner domain to its direct parent domain will be limited. In the vertical, 147 levels were used for all domains (80 levels within the lowest 2 km, 50 hPa model top; Figure 6c). Vertical grid spacing increases from ∼20 to 30 m within the lower 2 km and up to 1 km near the model top ( Figure 6c) A Rayleigh damping layer with an inverse time scale was used for the top 5 km to reduce errors from gravity waves reflected at the domain top [57], while sixth-order monotonic horizontal diffusion (coefficient value of 0.12) was applied to minimize spurious behavior at poorly resolved scales [58]. Other physical parametrization schemes used were the RRTMG scheme [59] for radiative transfer, the MM5 surface layer scheme [60], and the Noah land surface model [61].
For all simulations, the topography was based on the Geospatial Information Authority of Japan (GSI) 50 m digital elevation model (DEM). Simulations were carried out using the model's default terrain following coordinates, with a minimal amount of smoothing to ensure numerical stability over steep slopes. Even after smoothing, the DEM is able to represent the topographical features, such as the ridges of the west side of the volcano, well (Figure 6b; orange lines). As complex terrain is usually the reason for computational instabilities, which necessitate small time steps, a sensitivity simulation was also carried out for the June case using the default Global 30" elevation (GTOPO30) DEM (Figure 6b; gray lines). After this sensitivity run, we have set the time step for D4.
Because of the complex features of the model topography generated by the highresolution DEM data, the time steps of the temporal integration from D1 to D4 are set to 18, 6, 0.1, and 0.033 s, respectively.
The simulations were initialized using the ERA5 reanalysis data ( ∆x ∼31 km, 137 vertical levels). The outermost domain boundaries were specified using ERA5 hourly data. The simulations for the PBL domains for both cases started at 1200 UTC (2100 JST) of the previous day and were run for a total of 12 h to allow for the model spin-up. Simulations for the LES domains were initialized a later point to minimize computational time: six hours in the simulation for D3 (three hours of spin-up time) and eight hours for D4 (1 h spin-up time to allow for the main flow to travel from the lateral boundaries past the main study area). Data were output every minute.

Particle Trajectory Modeling
In order to quantify the magnitude of the impact the detailed airflow over the complex volcanic terrain has on the fall speed of particles near the volcano, idealized particle fall calculations were carried out over the 3 h average wind fields (0600-0900JST). Initially, 100 particles were randomly introduced in the domain using an ellipsoid distribution.
Idealized particles trajectories were calculated here. Calculations were carried out for two vertical velocity fields: (i) W 0 = −U T , where −U T is the terminal velocity of the particle, representing the case of particles falling unperturbed, and (ii) W Pert = w − U T , where w (the WRF vertical velocity output) was added to represent the influence of detailed topography. In both cases, calculations were carried out using the model output for the horizontal wind fields. Five terminal velocities were chosen: 8, 4, 2, 1, and 0.3 m s -1 , to represent typical sizes and types of ash particles based on the formula [5], corresponding to 1, 0.5, 0.25, 0.125, and 0.063 mm for single particles (density ρ A = 2500 kg m -3 ) and 2, 1, 0.5, 0.25, and 0.125 mm for an aggregate (ρ A = 750 kg m -3 ).
It is noted that the effects of heating of the atmosphere from volcanic eruption are not incorporated in the present trajectory modeling. How to formulate the effects of heating due to volcanic eruption in the trajectory modeling will be an interesting topic; however, this study does not take into account this heating effect, because the central focus here is placed on the influences of the complex topography on airflows and particle transport and sedimentation.

Stationary Flow Characteristics
Output from the innermost WRF model domain (D4) is shown as 3-h averages in Figure 7. In the June simulations (Figure 7a,b) flow is from the east (i.e., right side) up until the flow reversal point at 2 km. Mountain waves are triggered over the ridge to the east of the domain and in the lee side of the volcano, but are damped with distance due to the effect of the boundary layer acting as a sponge layer [62]. The location of the upwards component of the wave on the volcano lee matches the lidar observations (shown as black wind vectors). Gravity wave activity is also completely damped above the flow reversal point. A layer of flow reversal can be seen between the eastern ridge and the volcano up to 250 m from the surface, and a lee-side vortex can be seen immediately to the lee of the volcano. Finally, two isolated areas of flow reversal can be seen between the gravity waves and the westerly flow layer, from 1.6-2 km. Simulations using different DEMs are qualitatively similar. Gravity wave activity is also dominant in the October simulation (Figure 7c). Over the eastern ridge, gravity waves are triggered over the inversion layer, with two rotors forming between the waves and the inversion (seen as areas of flow reversal). Over the volcano, the flow response is a hydraulic jump. Lidar observations match with the position of the upwards component of the wave and show an increased magnitude over the June case.
The wind direction shear and multiple topographic maxima combination leads to a complicated surface flow for both days, as visualized by streamlines. Despite significant differences in the vertical structure above ∼1 km, common elements can be seen for both case study days (Figure 7d-f). Inside the Aira caldera, the general flow is to the southwest, interrupted by the volcano causing flow divergence. A stagnation point can be seen for both days in the west side of the volcano (at approximately x = 0 km, y = 3 km) and an extended area of lee-vortices can be identified in the south-western side of the volcano. In an idealized setting, flow stagnation points are expected to be along the same axis as the lee vortices, but this is not the case here, owing to the wind direction shear. As with the vertical cross-section, using a 30" DEM results in a qualitatively similar flow (Figure 7d,e).
The resolved flow is also in agreement with most surface observations that show a northeasterly flow. The location of the lee vortex is captured well, as seen by the station in the southwest of the volcano, which has a very small wind speed reading. Subtle changes in the wind direction owing to flow splitting are seen, although in the June case flow stagnation is overestimated in the western side of the volcano (∼x = 5 km, y = 1 km). The single worst-matching point is located over the southern part of Kagoshima (∼x = 10 km, y = −2 km), which recorded an almost constant southeasterly wind, in contrast to the northeasterly wind recorded by the other stations in the city.
At 600~m asl, the flow is less complicated, with southeasterly winds in the June simulations (Figure 7g,h) with a more pronounced easterly component in the October case (Figure 7i). Streamlines show significant diversion as the flow meets the volcano. On the lee side, the gap created by the diverted flow is filled with the downslope wind seen in Figure 7a-c. A comparison against the lidar data at this height shows that the average wind direction is captured well in the lee of the volcano.

Lee-Side Turbulence Characteristics
The position of the Doppler lidar in the lee of the volcano provides an opportunity to compare the way mountain waves were modeled and the related turbulence characteristics. As the lidar was only able to consistently provide observations over a thin layer, results presented are averages between 500-700 m asl (centered at 600~m asl). Figure 8 shows average values for the three wind components around the lidar side between 0600-0830 JST. Areas of significant absolute vorticity are also shown. A narrow wake region can be identified in the area between changing signs of absolute vorticity.
The wake region can also be identified as a local minimum in the u component. In the case of the two 50 m DEM simulations (Figure 8a,c), the wake is generated by the lee-side ridges via a streak: vorticity is generated over the ridge as small changes in the wind direction lead to a downslope flow to either the left or the right side of the ridge. In the 30" DEM simulations, the wake region is isolated (Figure 8b) and vorticity is generated due to the flow around the volcano (e.g., as in [63]). In the case of the June simulations, a series of weak vortices can be seen as changes between positive and negative values of both the v and w components (Figure 8d,g). The resolved wave pattern as seen through the averaged vertical velocity values (Figure 8g-i) closely mirrors the shape of the topography depending on the incident angle.
In the June (50 m DEM) simulation, a strict comparison against the lidar observations shows an overestimation of the u and v components (Figure 8a,d) and a correct representation of the w component (Figure 8g). On the other hand, most values are represented appropriately in the October case, with only a small underestimation of the v component (Figure 8f). In the case of less intricate topography (e.g., the volcano in the 30" DEM simulation), the lee side patterns of all wind components are considerably smoother (and provide a worse match between lidar observations). As such, considering the proximity of the lidar location to the lee-side ridges and their connection to the low velocity streaks, it can be assumed that such errors in the spatial representation are tied to the complexity of the terrain and are inevitable to a certain degree. The wake region is also associated with increased resolved turbulence kinetic energy (k) and low-period (high-frequency) oscillations in the vertical ( Figure 9). Specifically, maxima in the k distribution can be identified along areas of absolute vorticity sign change (Figure 9a-c). In the two 50 m DEM simulations, the area of increased k intersects the volcano, while in the 30" case it is disconnected from the topography. The June simulations show extended areas of k > 0.4 m 2 s -2 , while in the October case, high k values are only resolved across a narrow area in the lee of the volcano. Turbulence kinetic energy calculated based on the Doppler lidar measurements was 1.6 m 2 s -2 for the June case and 2.3 m 2 s -2 for the October case. The maximum resolved k in the model in both 50 m DEM cases is 1.5 m 2 s -2 , showing that a significant part of the observed wave-related turbulence was resolved. In the October case, the area of increased k overlaps with the Doppler lidar site; however, in the June case, the maximum values are located ∼1 km away from the site. This mirrors the spatial error seen in Figure 8.  (Figure 9d-f). However, areas of increased k tend to be associated with low-period oscillations, similar to those seen in the Doppler lidar signal (Figure 8). In the June simulations, areas of low period can be identified scattered across the volcano lee-side, while in the October case, they closely mirror the area of increased k. In both cases, observed maxima in the w distribution (10 and 3 min for the June and October cases, respectively) are resolved in the model, over (October) or near (June, 50 m) the Doppler lidar site.
As discussed in previous studies [64][65][66], the effects of turbulence on flow and dispersion are important in high-resolution simulations. In the present study, we intend to reproduce flow fluctuations generated by complex topography with the use of highresolution DEM data. From the results shown in Figure 9, we are able to capture the turbulent flows observed by Doppler lidar in a certain degree. Because this study does not employ turbulence generation method at the lateral boundaries like those used in a turbulence recycling approach [34], fully turbulent flows such as seen in flows over rough surfaces would not develop. Complete LES modeling with a proper turbulent generation method at the lateral boundaries may be necessary to represent fully turbulent flows.

Sensitivity to Resolution
The WRF model is able to reproduce the vertical structure of the atmosphere reasonably well, as seen by a comparison against atmospheric sounding observations ( Figure 10). The associated root mean square errors (RMSE) for all domains are small and generally decrease with an increase in resolution. At the innermost domain, the final RMSE values for wind speed and direction range between 1.2-2.4 m s -1 and 11-13 • , respectively (Figure 10a,b). Although the error in wind direction decreases with increased resolution for the October case, this is not true for the June case, where RMSE values oscillate from 11-13 • . WRF output also shows significant improvement in the representation of dynamic stability (shown by Richardson number, Ri; Figure 10c). In the June case, increased resolution leads to a small decrease in fidelity in relative humidity (1% increase over the outermost domains). Comparing the two June simulations shows a very small decrease in fidelity when using the 30" DEM (star markers). Note that calculations were carried out from 0-10 km and as averages over a 31 km area at a single time (0900JST). A comparison between average wind speed and wind direction values for each domain is also shown in detail in Figure 11. Increase in wind speed fidelity with increased resolution is mainly caused by the better representation of very low wind speed values, i.e., less than 2.5 m s -1 (Figure 11a-d). In the June case, as most of the wind speed observations are small and as such resolved better in the innermost domain. On the other hand, the October case features a wider range of wind speed values, between 1-5 m s -1 . Increased model resolution in this case leads to a better match over the lower end of the range but underestimation at some points in the higher end (i.e., between 2.5-5 m s -1 ). Changes in wind direction do not exhibit a similar systematic improvement (Figure 10e-h). In all cases, wind direction away from the volcano is associated with small errors, while the stations over the volcano present a more significant challenge, especially in the June case that features overall lower wind speeds. As the surface flow has the largest amount of parameterized turbulence, results here point towards the existence of a minimum error with the current modeling approach.
Note that RMSE values shown for the AEROS stations (i.e., the instantaneous value output stations) are the minimum RMSE for that point in the domain at the time of the observations ±7 min. The time window was chosen to allow for a fair representation of the inherent turbulence that is included in the observations and was decided as the average time period for which the distribution of the perturbation part of the horizontal wind component differed less than 10% when compared to a 1 h average (i.e., the minimum time required to get a representative sample of randomness in the model output).

Particle Fall Trajectory Analysis
Overall, considering the complexity of the topography and vertical structure the model performs well even at coarse resolutions. Results show a consistent increase in fidelity and the range of w was shown to significantly increase for the innermost domain. Although a comparison with Doppler lidar measurements revealed spatial inconsistencies, the resolved values were seen to be realistic for all three wind components. As such, results from D4 are used here to provide a basis for the idealized particle fall trajectory calculations of small to large particles (U T = 0.3, 1, 2, 4, and 8 m s -1 ). As a first step, 3-h average wind field values are used to provide an order of magnitude evaluation. More precise calculations necessitate the use of a specialized transport model, further complicating the study. As such, simulating ash dispersal for the events will be the point of a future effort.
Example trajectories for the W 0 (no vertical velocity component) and W Pert (WRF vertical velocity added) are shown for an injection height of 500 m and three particle categories in Figure 13. In all cases, including the vertical velocity component leads to visibly different trajectories. In the W Pert calculations, large particles (U T = 8 and 2 m s -1 ; Figure 13a,b,d,e) are captured by the downslope wind, which holds the trajectories close to the volcano's slope (as seen in Figure 2b) and enforces deposition closer to the release location. In the intermediate case (U T = 2 m s -1 ; Figure 13b,e), inclusion of the vertical velocity component leads to markedly different trajectories, with the deposition of particles occurring over the entire lee side of the volcano. The effect of the wind shear is also visible in all calculations, in a similar manner to Figure 2a. Finally, in the case of small particles, gravity wave activity mainly results in the vertical displacement of particle trajectories. Between the two days (first and second rows in Figure 13), results are similar; however, the differences between the W 0 and W Pert calculations are emphasized in the October case due to the larger horizontal wind component. The impact of the gravity waves on the fall speed is shown via the ratio between the average fall speed to the terminal velocity (i.e., the unperturbed fall velocity), with values away from 1 showing a large deviation ( Figure 14). In this sense, the impact is most notable for intermediate particles (U T = 1-4 m s-1; medium-sized ash particles). A comparison between the June and October simulations (Figure 14a,c) shows that the effect is emphasized in the case of more energetic flow conditions (i.e., the October case), while a comparison between the 50 m and 30" DEM simulations (Figure 14a,b) shows that the smoother flow simulated using the 30" DEM leads to a decreased impact on the fall speed. Finally, the injection height of the particles can also have a direct impact: in all cases, injecting particles at a 1000 m agl leads to ratios closer to 1. A secondary characteristic is the range exhibited due to the random nature of the particle placement, which increases significantly for smaller terminal velocities (i.e., smaller particle sizes). Specifically in the case of U T in the June simulation (Figure 14a), although the ratio is close to 1, the standard deviation shows ratios constrained within a factor of 2, with values less than 1 showing particles that fall slower than their terminal velocity. This is due to the vertical displacement of the particle trajectories that was previously seen when discussing Figure 13c,f.
Finally, the impact of detailed airflow can be seen as a change in the deposition distance of the particles ( Figure 15). As the perturbed fall velocities are, on average, higher than the terminal velocities, the deposition distance is generally smaller in the case of the W Pert calculations. For an injection height of 1000 m (Figure 15a-  In the case of an injection height of 500 m (Figure 15d-f), the difference in deposition distance becomes more pronounced for lighter particles, i.e., U T < 2 m s -1 , with a maximum difference of 3 km in the June case (U T = 1 m s -1 ) and up to 4 km in the October case (U T = 2 m s -1 ). Furthermore, in both June and October cases, 2-9% of the particles were forcedly deposited due to downslope winds; at U T = 0.3 m s -1 for the June case and U T = 1 m s -1 for the October case. In both cases, the average deposition distance is similar to that of the immediately larger particle size category (i.e., 0.3, 1 m s -1 for the June case and 1, 2 m s -1 for the October case).

Discussion and Conclusions
The WRF model was used in a multi-scale setting in order to conduct high-resolution simulations for two cases with a similar, complex vertical structure in the atmosphere. The Sakurajima volcano, Japan, was used as the study location due to the observational network around it and the need of detailed simulation data in order to understand the dynamics of transport of pollutants from the volcano. Four domains were used, with the outer domain initialized ERA5 reanalysis data. The two outermost domains employed a PBL parameterization, while the innermost two employed an LES scheme. The atmospheric motions resolved were seen to have two major components, a quasi-steady state response dominated by gravity waves, overlaid with localized turbulent perturbations in the wake regions of the various topographic maxima.
The overall steady-state flow pattern was resolved well in all domains, agreeing both with theoretical expectations (e.g., flow splitting, mountain wave generation, lee-side wake) as well as surface, near-surface, and synoptic data. Although the majority of surface wind observations were represented well in the model, misrepresentation of wind direction was also seen over some of the observation locations. In the most severe cases, a 120-140 • error in wind direction persisted even at the innermost domain. Misrepresentation of surface values of wind direction is a known issue with WRF: it was shown [67] that, especially in the case of complex topography, the associated error of wind direction correlates with wind speed-low wind speeds leading to large wind direction errors, which was also partially true for simulations shown here. Results here show that issues raised for simulations over coarser resolutions can be maintained even after the resolution allows for a faithful representation of the topography, highlighting limitations of the approach used.
Turbulent fluctuations were studied in detail for the lee-side of the volcano and it was seen that a realistic orographic flow response was resolved in the innermost (LES) domains. On average, the model was able to reproduce the observed characteristics of the flow response, both at the surface (flow splitting and surface rotors) and aloft (wake region). A comparison against Doppler lidar observations showed that a qualitative match was generally achieved in the simulations (i.e., location and magnitude of the various orographic phenomena); however, issues were noted concerning the spatial fidelity of the representation. The complexity in the topography of the western side of the volcano (lee-side in the cases studied here) is suspected to be a main cause of the mismatch.
Turbulence aloft was increased over a wake region and tied with gravity waves (a well-known connection); via changes in the absolute vorticity of the flow [68,69]. Turbulence generation was seen to occur as a combination of the large-scale flow around the volcano (resolved both in simulations using the detailed 50 m DEM as well as the smooth topography of the 30" DEM), overlaid with a region of increased turbulence caused by the small-scale changes in the flow over the ridges on the lee-side of the volcano (resolved on when using the 50 m DEM). Due to the sensitivity of the orographic response to small changes in the wind speed and wind direction [70], it is assumed that over highly complex topography, a perfect spatial match would be difficult to achieve without data assimilation. Simulations using the 30" DEM were shown to produce a smoother response on the leeside. As such, it can be assumed that a spatial match might be more easily achievable in simulations over a location with smoother topography.
Results from the 30" DEM simulation, although qualitatively very similar to those of the high-resolution (50 m) DEM, showed that the correct representation of topography can have a large impact on fidelity. Simulations using the 50 m DEM showed a consistent increase in fidelity in the representation of surface wind direction (by 6.5 • ), using the 30" DEM led to a deteriorating match, with an increase in the associated error for wind direction (by 10.7 • ). Smoothing techniques are often employed in computational modeling over complex terrain, but results here show that taken to a relative extreme (as in the case presented here) smoothing can have a deteriorating effect on the output of high-resolution simulations. Similar results were noted in a model inter-comparison study [71], which studied rainfall and wind patterns over Kyushu in typhoon conditions and concluded that the way topography is resolved in the models had a significant effect in the results. The importance of using high-resolution DEM data in high-resolution atmospheric simulations was also emphasized in a heavy rainfall case [72] and in a strong wind case [73].
The proper representation of the detailed airflow is known to play a vital role in the proper simulations of pollutant transport [18,45,73,74]. Under stable conditions (common for Sakurajima volcano; [45]) and for particles that are either introduced close to the surface or, in the case of large eruptions, have fallen from a large height and are captured by orographic flows generated by distal mountains (for example as in the case discussed in [75]), the impact of localized airflows was noticeable even for heavy particles (∼20% change in speed and up to 1 km difference in deposition distance). For particles that were not directly deposited, gravity wave activity was shown to force vertical displacement. As expected, the impact was seen to lessen with an increased injection height and can be expected to be negligible for particle transport a few kilometers above the surface. The calculation of particle settling velocity is a critical process in VATD models [12] and results here illustrate that modifying the calculated velocity over topographic maxima would be an easily implementable approach to parameterize the impact of orographic effects and improve accuracy both in VATD models, as well as backward trajectory calculations that are often performed on the field [75].
The main reason why this study used of an off-line model in the trajectory calculations is that in a volcanology viewpoint we would like to conduct a large number of sensitivity tests and follow a probabilistic approach due to uncertainties in the eruption source conditions. A fully coupled atmosphere-chemistry-aerosol model such as the WRF-chem performs favorably for a single deterministic simulation (and as emissions are handled internally it makes full use of WRF calculations on subgrid effects), but repeating the meteorological calculations after that can be significantly expensive in most volcanological settings (especially in the case of LES). Thus, the present off-line modeling framework should be useful in examining possible flow and dispersion patterns with some assumed eruptions conditions.
Overall, results further highlight the difficulties of accurately simulating volcanic ash deposition in the proximal region. Even using the modeling presented here, spatial fidelity was shown to be limited within a factor of O(100 m) due to the sensitivity of the resulting orographic effects on the inherent errors in the initialization data [70] as well as due to unresolved turbulence. In the case of VATD modeling, meteorological data are used along with estimated eruption source parameters (e.g., erupted mass, volcanic plume height, duration, ejected particle grain size) in order to estimate both the plume characteristics (e.g., shape, radius [76], as well as the transport of volcanic ash [77]). This leads to such applications being especially sensitive to systematic errors in the meteorological data as it affects both the initial condition (i.e., plume shape) as well as the transport. In order to counteract this effect in cases of high-resolution secondary modeling applications, an ensemble method could be used with slightly different initialization times to account for the randomness in the turbulence, or in cases where this is inapplicable, an empirical increase in horizontal and vertical diffusion. Recognizing these issues will contribute to deepen the understanding of atmospheric transport of volcanic ash with the use of a fully coupled atmosphere-chemistry-aerosol model such as the WRF-chem modeling system [18,78].