Use of Thermal Image Velocimetry to Measure a Dust-Devil-Like Vortex Within a Sports Ground in a Residential Area

We investigate the characteristics of a dust-devil-like vortex (DDLV) observed using thermal image velocimetry (TIV) at a sports ground in Tokyo. Thermal image velocimetry provides unique observations of the two-dimensional velocity distribution for a DDLV with high spatio–temporal resolution (i.e., tens of cm s−1) near the ground. Two DDLVs were detected, one each in summer and winter, and the quantitative features of the larger, stronger DDLV in the winter are examined. The size and strength of the detected DDLV, which are quantified using TIV, are within the ranges reported in past observations and numerical simulations of dust devils. The vortex appears at the boundary of a cold-air current near a 55-m building wall, and persists for more than 3 min.


Introduction
Dust devils are thermally-induced vertical vortices comprising sediments that are picked up from the ground. They are usually on the order of 10 m in diameter, and 10 m s −1 in tangential velocity (Stull 1988), and persist for a few minutes. Murphy et al. (2016) summarized ground observations of dust devils with sizes ranging mostly from 1 to 100 m in diameter and lifetimes from seconds to minutes. Dust devils account for a considerable amount of sediment pickup in deserts and bare-soil areas, which directly contributes to global-scale transport of dust sediments (e.g., Uno et al. 2005). It is also likely that dust devils control the opacity of the Martian atmosphere (Balme and Greeley 2006). Fujiwara et al. (2011) observed many dust-devil-like vortices (DDLVs) in an urban area using Doppler lidar. Dust-devil-like vortices are vertical, swirling vortices that are as strong and large as dust devils, but contain no dust due to the absence of a dust source on the ground B Atsushi Inagaki inagaki.a.ab@m.titech.ac.jp 1 Department of Transdisciplinary Science and Engineering, Tokyo Institute of Technology, Tokyo 152-8552, Japan (e.g., Kanak 2005;Fujiwara et al. 2011). Raasch and Franke (2011) reported many DDLVs in a convective mixed layer. Dust-devil-like vortices are considered a ubiquitous feature of convective mixed layers (e.g., Ito et al. 2010Ito et al. , 2013, and have a significant influence on vertical momentum and heat transport near the ground (Kaimal and Businger 1970).
Some researchers have considered the initiation of dust devils and DDLVs to be associated with vertical airflow shear near the ground, due to opposing warm and cold airflow lifting a horizontal vortex tube near the surface (e.g., Renno et al. 2004). They are associated with multi-directional convergent flows on the ground below convective cells (Raasch and Franke 2011). The surface-related process is more complex in cities due to buildings, which act as obstacles making the creation of coherent horizontal vortex tubes and/or horizontal vortex sheets difficult. Another initiation mechanism over very rough surfaces is attributed to flow separation at the edges of obstacles (e.g., Williams 1948). Coceal et al. (2007) reported that roll-up of a strong shear layer behind cubical obstacles initiated a pair of hairpin vortices, which further evolved into very large streak-like structures within the logarithmic layer. Fujiwara et al. (2011) reported DDLVs located 150 m above the ground in an urban area, although they concluded that the location of DDLVs at that height is not directly related to urban geometry. Remote sensing can be used to detect the velocity field within and around dust devils and DDLVs. Their two-dimensional velocity distribution has been obtained via ground-based remote sensing at a fine resolution. Specifically, Bluestein et al. (2004) used Doppler radar to measure the velocity distribution within dust devils 5-10 m above the ground, observing their variability over the course of tens of seconds. Fujiwara et al. (2011) used Doppler lidar to detect DDLVs and their surrounding flow distributions. However, Doppler radar can detect the velocity field only within dust devils, and not in their surroundings, where there is no aerial dust. In addition, there is also a limitation in the measurement height for Doppler radar due to significant reflection from the ground. Although Doppler lidar on the ground surface is potentially able to detect near-surface wind, it is still difficult in locations with very rough topographies, such as built-up areas. Both single-Doppler radar and lidar can measure the radial velocity distribution, but cannot measure two-dimensional velocity vectors.
Regarding flow visualization, Vogt (2008) was able to detect flow distribution over a wheat field at a very fine scale of~100 m based on time-sequential thermography, which can extract the footprint of unsteady near-surface temperature structures (e.g., Paw et al., 1995;Christen et al. 2011). Inagaki et al. (2013) combined this method with particle image velocimetry to measure the near-surface velocity distribution. This thermal image velocimetry (TIV) technique has been used to measure two-dimensional velocity distribution 5 cm away from a vertical building wall.
In this study, we analyzed a DDLV, which was observed using TIV within a sports ground surrounded by buildings in Tokyo, Japan. The TIV measured the two-dimensional distribution of the velocity field within and around the DDLV with high spatio-temporal resolution (< 0.3 m in 1 s). This rich dataset was used to quantitatively evaluate the properties of the DDLV, as well as its temporal evolution.
It is mentioned that the DDLVs were observed next to a building taller than 50 m, and therefore might be related to the building wake. However, the DDLVs that we observed were exceptionally larger and had longer lifetimes, and were also much more rarely observed compared to the other vertical swirls. Therefore, we focus on their characteristics as rare but significant events possibly occurring in urban airflow.

Thermal Image Velocimetry
Thermal image velocimetry detects near-surface wind velocities by tracking the surface temperature pattern. The technique assumes that the surface temperature pattern represents the footprint of a coherent flow structure moving with the wind. Inagaki et al. (2013) provide a detailed description of TIV.
Surfaces that are spatially heterogeneous in terms of emissivity, albedo, heat capacity, and geometry produce a heterogeneous surface temperature pattern, without heterogeneity in the wind structure. Vogt (2008) used short-term temporal averages to separate the stationary pattern attributed to heterogeneous surface properties from the variable pattern attributed to the coherent structure of winds, as follows where T s is surface brightness temperature, and the overbar and prime symbols denote the short-term temporal average and deviation from the average, respectively. The deviation component has been shown to accurately describe the spatial wind pattern (e.g., Vogt 2008;Christen et al. 2011;Garai and Kleissl 2011;Garai et al. 2013;Inagaki et al. 2013;Morrison et al. 2017).
The temporal filter size is implicitly determined by the size of eddy motion being tracked by the TIV. Therefore, it depends on experimental and analytical parameters such as image resolution, sampling frequency, and the interrogation window size for the image correlation analysis. These parameters limit the detectable size of the turbulent footprint. It is also dependent on the thermal admittance of the target materials (e.g., Christen et al. 2012). An object with higher thermal admittance requires a longer period of heat accumulation or dissipation to change its surface temperature enough to be detectable by the thermal infrared (TIR) sensors. Such an accumulation period is accomplished during the passage of an eddy structure with a longer lifetime whose signature can be effectively extracted using a longer filter. The specific settings for the analytical parameters were determined based on the study by Inagaki et al. (2013).

Measurements
Surface brightness temperature was measured within a sports ground at the Tokyo Institute of Technology, Japan (35°36 17.6 N, 139°40 55.6 E). The sports ground is surrounded by residential houses to the west and south sides, and relatively tall university campus buildings (approximately 20-50 m in height) to the north and east sides. The ground is covered by 30-mm-thick artificial turf.
We used a TIR camera (SC5200; FLIR Systems, Wilsonville, Oregon, USA) to measure the surface brightness temperature pattern. This camera captures thermal infrared radiation for wavelengths of 2.5-5.1 µm at a frame rate of 100 Hz. This device system can detect temperature fluctuations (sensitivity of 0.002 K) with low background noise, which is suitable for application to TIV. The same type of device was used by Morrison et al. (2017) for thermal footprint detection in a desert area. The camera was installed at the top of the W8 building (hereafter, W8) at the Tokyo Institute of Technology located to the east of the sports ground, at a height of 55 m above the ground. Figure 1 shows the location of the camera and its  field-of-view (FOV). The measurement variables are defined following the Cartesian coordinate system, in which x and y axes are set as shown in Fig. 1. Five sonic anemometers (DA600 TR90-AH; KAIJO, Tokyo, Japan) were installed at different heights within the FOV of the TIR camera, namely, at 0.3, 0.6, 0.9, 1.2, and 1.6 m, all aligned vertically. Measurements from the anemometers were used to validate the velocity measurements from the TIV, and to examine local atmospheric conditions when DDLVs occurred. In addition, two thermocouples were installed within the canopy of the artificial turf near the sonic anemometers. Data from the sonic anemometers and thermocouples were recorded using a data logger (LR8400; HIOKI E.E., Nagano, Japan) at a 50-Hz sampling rate, all in sync. We conducted measurements for about 1 h at approximately noon on 10 days (1 day in February and 9 days in September in 2013). These days had clear skies, except for occasional cumulus clouds partially shadowing the ground.
Additional meteorological datasets to help describe the general weather conditions were obtained at the top of W8, and at the M1 building (hereafter M1), located about 400 m apart in the northern part of the same university campus. Single weather stations (WXT, Vaisala Oyj, Helsinki, Finland) and three-dimensional sonic anemometers/thermometers (Model 81000; R. M. Young Co., Traverse City, MI, U.S.) were installed at the W8 and M1 sites, respectively. We also referred to meteorological datasets obtained at Otemachi (hereafter, OT) by the AMeDAS observatory, which is operated by the Japan Meteorological Agency. The OT site is over 11 km away in the north-east direction. Table 1 summarizes all instruments, locations, and the variables measured at each site.

Weather Conditions
The observations were conducted on clear-sky days when sunshine duration was mostly 100%. During the 10 observation days (i.e., 20 February, and 2, 3, 9, 12, 13, 17, 19, 27, and 30 September 2013), DDLVs were observed twice at 1200 LT (local time = UTC + 9 h) in the Japanese winter and summer. We focus on the DDLV that was observed at approximately 1200 LT on 20 February 2013 (see Online Resource 1). This day had clear skies with partial cumulus-cloud cover. The general weather conditions at this time are described by measurements from the four observation sites (i.e., W8, M1, OT, and the sports ground). We focused on the DDLV in winter because it was larger and stronger than the one in summer. In addition, the summer DDLV was in the video frame for less time and the measurements were partially obstructed. However, the general properties of the DDLVs in winter and summer did not differ significantly, as shown in the Appendix. Figure 2 shows the diurnal variation of mean air temperature, wind speed, and direction for our case study. The measurements at M1, W8, and the sports ground are averaged every 1 min, while those at the OT site are averaged every 10 min. Although the measurement points are physically separated in space by tens of metres (up to > 10 km), the diurnal patterns of each variable correspond well to each other. Therefore, these measurements represent the general regional weather conditions.
The environmental conditions when the DDLV was measured are described as follows. A DDLV was observed when: the air was unstably stratified with a temperature difference of more than 10 K between the near-surface air temperature at the top of the turf (T s ) and that above the rooftop of buildings (i.e., 55.1 m above the ground, W8) (Fig. 2a), and the wind speed was less than 1 m s −1 at both 1.6 m and 55.1 m from the ground (Fig. 2b). This wind speed is lower than that of Fujiwara et al. (2011). This period of low wind speed (i.e., 1100-1200 LT) is associated with the transition of regional flow from a northerly to a southerly wind (Fig. 2c). These conditions, especially the weak horizontal wind speed and unstable temperature stratification, favour dust devil occurrence. The atmospheric stability was also measured using the sonic anemometers within the sports ground, using the following where z is height from the ground, L is Obukhov length, g is acceleration due to gravity, w T is kinematic sensible heat flux, u * is friction velocity, and k is the von Kármán constant. The friction velocity was measured using a sonic anemometer at the site based on the eddycovariance momentum flux. The value of z/L was − 1.42 (dimensionless) at 1.6 m from the ground at the location and time of the observed DDLV, implying strongly convective conditions. The ratio of the convective velocity scale w * to the friction velocity is also used to determine necessary conditions for the occurrence of the dust devils (Rafkin et al. 2016). The convective velocity scale is defined as follows where z i is the atmospheric boundary-layer height. Although we have not directly measured z i , it was estimated as 309 m, which is the average of the direct measurements using Doppler lidar (Yagi et al. 2017) at the same site from 1 January 2013 to 11 February 2013. As a result, w * /u * = 8.8 for our winter observations, which satisfies the necessary condition of being greater than 5. The air temperature less than 10°C is much lower than that during the past observations of dust devils (e.g., Rafkin et al. 2016). One exception is reported by Fujiwara et al. (2011), who measured DDLVs when the air temperature was less than 10°C at 1.5 m above the ground, observed at the AMeDAS observatory in Sapporo. Since they clearly showed the relationship between DDLVs and thermal convection structures, the mechanism that sustains the swirling motion is considered to be similar to that for dust devils.

Thermal-Image-Velocimetry Analysis
We examined TIV flow and velocity outputs and applied them to our DDLV analysis. The surface brightness temperature pattern (Fig. 3a) is mainly determined by the colour of the artificial turf (Fig. 3d). After subtracting a short time average (5 s) from each pixel, a clear swirling structure was visible, containing elongated temperature "streams" centred at (0, − 30) m (Fig. 3b). By tracking the temperature fluctuations within an interrogation window size of 15 × 15 pixels among two sequential images separated by 0.1 s, we derived a vector distribution (Fig. 3c). Visual inspection of the processed images revealed that the swirling motion was observed for more than 3 min before the TIR camera stopped recording, the swirl had a diameter of about 30 m, and the swirl rotated in an anticlockwise direction.
When the mean flow was in the opposite direction, another DDLV was observed, this time rotating in a clockwise direction (see Online Resource 2). This implies that the direction of a DDLV is insensitive to the direction of the Earth's rotation, as noted in Murphy et al. (2016), and is more sensitive to the wind direction and the geometries of surrounding obstacles. This is further examined in Sect. 4.2. Figure 4 shows a time series of the advection velocity of the surface brightness temperature pattern (hereinafter TIV velocity), and the wind velocity as measured by a sonic anemometer at a height of 0.75 m above the ground during the period in which the DDLV was observed; they are averaged every 1 s. The velocity components U and V are parallel to the x and y axes shown in Fig. 1, respectively. The variation in wind speed and TIV velocity correspond well. Figure 5 shows the vertical distribution of mean wind speed, coefficient of correlation, and slope of the linear regression line between the wind and TIV velocity, u SA and u TIV respectively, at each measurement height. The bottom measurement value of u SA exceeded the value at the second measurement level (z = 0.3 m). This is probably due to the flow distortion by the other instruments or their supports, although the fluctuation signals were still considerably correlated with u TIV . The coefficient of correlation exceeded 0.8 at all measurement heights (Fig. 5b). The ratio of u SA to u TIV becomes unity at around 0.8 m above the ground ( u SA u TIV = 1, Fig. 5c). Therefore, the TIV velocity represents the wind velocity at 0.8 m from the ground. Figure 6 shows the time series of the ratio of u TIV to u SA for different heights. The mean values are steadily constant during this observation period, although they fluctuate in time, and the standard deviation is 0.25 for all measurement heights. A considerable deviation was This might be due to the passage of localized strong turbulence or the transient motion associated with the changing regional wind. We do not investigate it further in this study since it is well apart from the main vortex core.

Mean Structure of the Dust-Devil-like Vortex
To examine the mean structure of the observed DDLV, flow fields within and around the vortex were temporally averaged by tracking the vortex core, and spatially averaged in the tangential direction. The core of the DDLV was detected (1) subjectively by visual inspection of TIV images showing the velocity and distribution of brightness temperature, and (2) by finding the point of maximum correlation between the Rankine-vortex model and TIV velocity distribution, following the method of Fujiwara et al. (2011). Using method (1), we detected the vortex core for the period 935-1104 s after the observations started. Although the centre Fig. 6 Time series of the ratio of mean wind speed to advection speed for different heights was outside the camera's FOV (i.e., the period from 969 to 999 s), the position of the core was roughly estimated from the velocity distribution in the vortex sleeve. Method (2) requires the core to be inside the FOV of the camera, so it detected the core within a more limited time period of 1060-1099 s. As shown in Fig. 7, the track of the DDLV core detected by method (2) overlaps well with that of method (1). The DDLV appeared in the frame at 935 s when the centre of the DDLV was 20 m away from the wall of W8 to the east. It moved southward, turned in a clockwise direction for about 30 s (from 969 to 999 s), and finally moved northward. This path followed the prevailing wind, as shown in Fig. 2. The average speed of the DDLV was about 2.7 m s −1 , which is larger than the background wind speed at a height of 1.6 m above the ground (see Fig. 5). This is probably because the DDLV is taller than 1.6 m and is advected by the wind at a higher elevation. Figure 7 also shows that the DDLV initially appeared near the north edge of the W8 building, which is the largest building surrounding the measurement area. Ordinally, mean flow is separated from a building wall and generates vertical swirls in a wake motion. This is a possible mechanism for triggering a DDLV in a built-up area, as indicated by Williams (1948). The mean wind direction and the rotation direction of the DDLVs seen in winter (northerly flow, anticlockwise) and summer (southerly flow, clockwise) also support this mechanism. Normal building wakes cannot persist very long as DDLVs, especially under unstable stratification. Also, building wake phenomena are much more common than DDLVs. Therefore, the persistent mechanism of DDLVs observed in this area cannot be explained solely by building wakes. Regarding the initiation mechanism, the observation was made on a cold winter day, when the convective structures in the mixed layer were probably too weak to create strong convergence lines near the surface and roll up the vortex sheet (e.g., Renno et al. 2004;Raasch and Franke (2011). This is expected because most dust devils are observed on hot summer days over flat fields (Rafkin et al. 2016). Building wakes are an alternative initiation mechanism for winter DDLVs in urban areas, although this has to be examined with more observations. Figure 8 shows the radial distribution of mean tangential and radial velocity components, vertical vorticity, and temperature fluctuation averaged over 40 s around the vortex core, as detected by method (2). The peak tangential velocity component occurred 3.6 m away from the centre of the swirl, and was 4.2 m s −1 . This is slower than the criterion of Sinclair (1973) who suggested it should be on the order of 10 m s −1 , although it is similar to that of Kaimal and Businger (1970), of 5 m s −1 . Similar values were also obtained in a highresolution numerical simulation by Raasch and Franke (2011), who attributed this to the ensemble average decreasing in magnitude, rather than being an instantaneous value. The other variables were also qualitatively similar to those of Raasch and Franke (2011), especially the radial distribution of the radial velocity component. The peak vertical vorticity was 1.8 s −1 , and the core radius of the swirl was about 12 m, judging by the vertical vorticity profile and the distance between the centre and the point where the vertical vorticity became zero. Some research indicates that the radial profile of the outer sleeves of dust devils follows the modified Rankine-vortex model, described as follows where U θ and U θ m are the mean tangential velocity component and its maximum value in a vortex, r is a radius from the vortex centre, r m is the radius of the vortex core, and α is the decay factor controlling the decay of the vortex sleeve. The decay factor expresses the effect

Temporal Evolution of the Dust-Devil-like Vortex
The temporal evolution of the DDLV was visualized using two-dimensional streamlines ( Fig. 9), which are unique observations accomplished using TIV. Before the vortex core appeared inside the measurement area (Fig. 9a), a northerly wind penetrated into the predominantly westerly wind field. The frontal boundary of the northerly wind, which is identified by the densely-gathered streamlines, gradually progressed to the south (x = 0 to 10 m in Fig. 9b, c). There is significant horizontal shear along the front, which tends to generate vertical vorticities having the same rotation direction as the observed DDLV. After that, the whole frontal shape became rounded (Fig. 9d) and the DDLV appeared at the edge of the front (x, y) = (7, − 68) m (Fig. 9e). Finally, the DDLV turned to the north (Fig. 9f, see also Fig. 7 or Online Resource 1) almost simultaneously with the ambient background flow, which also turned to the north (Fig. 2c). Based on our observations in Fig. 9, the DDLV appeared on a convergence line of two different air masses moving in nearly orthogonal directions to each other, whose orientation likely determined the direction of rotation. The rounded frontal shape (Fig. 9d) might be attributed to the presence of tall buildings on the eastern side, which produced drag in the northerly wind that created anticlockwise vorticity. In urban areas, strong eddies are usually generated at building edges (e.g., Coceal et al. 2007).
The present two-dimensional-streamline analysis visually connects the frontal convergence line and the DDLV. The characteristics of the frontal convergence line were further examined using the ground-observation datasets. Figure 10  from the north, a decrease in air temperature of more than 3 K was observed. Similar cold currents have been detected using time-sequential thermography, by Vogt (2008) and Christen et al. (2011) in atmospheric surface layers. Therefore, these cold fronts are a ubiquitous feature of the atmospheric boundary layer, probably associated with the convergence line of the thermal convection seen under daytime calm-wind conditions. Generally, near-surface horizontal convergence can enhance vertical vorticity since it causes updrafts with vortex stretching in the axial direction (Rafkin et al. 2016).
The horizontal shear at the convergence line can also enhance the vertical swirls if they have the same rotation direction. The penetration of the mixed-layer structures to the ground is visible within urban and vegetation canopies, as shown by some large-eddy simulations (Inagaki et al. 2012;Patton et al. 2016). Our experiment only observed the DDLV immediately after its initialization and it is still not confirmed how the cold front relates to the initialization mechanism.

Conclusions
This study investigated a DDLV observed using TIV at a sports ground located in a residential area of Tokyo. The method was able to detect the DDLV from the two-dimensional velocity distribution near the ground surface, at a high spatiotemporal resolution. Based on the observations, the motions and vorticity were quantified as an ensemble structure. The vortex appeared at the front between two colliding air masses, which propagated in nearly orthogonal directions to each other and had a considerable temperature difference (~3 K). At the frontal boundary strong horizontal wind shear was generated and the DDLV also appeared at the edge of this frontal boundary. This DDLV persisted for more than 3 min, much longer than most vertical swirls with lifetimes of less than a few seconds.
The DDLV was observed on a cold winter day when the convective structures in the mixed layer might have been too weak to invoke DDLVs because most dust devils are observed on hot summer days over flat surfaces (Rafkin et al. 2016). Williams (1948) proposed an alternative initiation mechanism with DDLVs being caused by the wakes behind surface obstacles. Since there are few studies of DDLVs in urban areas, more measurements or numerical simulations are necessary to clarify the major initiation process of DDLVs and the relationship between DDLVs and building structures in such areas.
This study showed the applicability of TIV to measuring the wind field around a DDLV and provided unique spatiotemporal measurements of a DDLV. This method is promising for measuring near-surface airflow on other planets, because it is not necessary to install any devices at the ground surface; airflow can be detected using satellite TIR imaging only.
Acknowledgements This work was supported by JSPS KAKENHI Grant Number 20H02253, and the Tokyo Tech. Young Investigator Engineering Award, Tokyo Institute of Technology. We would like to thank Mr. Mutaz Arif, Ms. Yiqun Zhang, and Mr. Eiji Iwatsuka for their technical support with the experiment.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix: Dust-Devil-Like Vortex Observed on 12 September 2013
Another DDLV was observed at noon on 12 September 2013 is briefly introduced here. The wind speed = 0.36 m s −1 from the south-east, and the 30-min averaged value of z/L = − 1.1 at 1.6 m above the ground. Figure 11 shows (a) the filtered surface brightness temperature and (b) velocity vectors, which were measured using TIV with the same analytical parameters used in the winter case. The figure axes are oriented the same as in Fig. 3. A large swirl rotating in a clockwise direction was seen at around (x, y) = (−48, 25) m. It passed through the FOV from the south-west to the north-east for 1 min at a speed of about 1.7 m s −1 . The temporal evolution of this DDLV is shown in Online Resource 2.
The velocity and temperature distributions surrounding this vertical swirl are averaged in time and radius by centring the vortex core (Fig. 12). The core diameter is about 6-10 m, which is a similar size to that seen in the winter. The mean structure of this vortex is similar to that observed in winter, although the tangential wind velocity component and vertical vorticity are smaller.