Observational analysis of a waterspout sighted in Macao, China, on 1 June 2021

Waterspouts are a relatively rare, yet potentially impactful, meteorological phenomenon. They have been observed and documented around the world for at least a century (e.g. Gay, 1921; Hurd, 1928; Golden, 1974a,b; Miglietta, 2019). In many cases, waterspouts (or their over-land equivalents) form from non-supercell storms (Wakimoto and Wilson, 1989; Wakimoto and Lew, 1993; Van Den Broeke and Van Den Broeke, 2015; Miglietta, 2019). However, a ‘waterspout’ is defined as a tornado occurring over water (Glickman, 2000), including supercell-spawned cases. Waterspouts occurring in the Pearl River Estuary of southeast China have been studied using radar, LiDAR, surface observations and numerical modelling. Stronger waterspouts have been associated with supercell storms, but may or may not have been associated with the mesocyclones of those storms. Chan et al. (2020) analysed data from a single radar and compared it with numerical models to study the evolution and structure of a waterspout-producing storm exhibiting a mesocyclone. Hon et al. (2021) studied the ability of numerical simulations to predict waterspouts by comparing simulated and observed radar data from a storm producing a weak waterspout near Hong Kong airport. Chan et al. (2021) employed dual-Doppler and Ground-Based Velocity Track Display (GBVTD) techniques to LiDAR data, and included surface observations, to analyse the behaviour of very weak anticyclonic waterspout and similar vortices near to, and impacting, Hong Kong International Airport, evaluating the role of horizontal shearing instability in vortex formation,

Waterspouts occurring in the Pearl River Estuary of southeast China have been studied using radar, LiDAR, surface observations and numerical modelling. Stronger waterspouts have been associated with supercell storms, but may or may not have been associated with the mesocyclones of those storms. Chan et al. (2020) analysed data from a single radar and compared it with numerical models to study the evolution and structure of a waterspout-producing storm exhibiting a mesocyclone. Hon et al. (2021) studied the ability of numerical simulations to predict waterspouts by comparing simulated and observed radar data from a storm producing a weak waterspout near Hong Kong airport. Chan et al. (2021) employed dual-Doppler and Ground-Based Velocity Track Display (GBVTD) techniques to LiDAR data, and included surface observations, to analyse the behaviour of very weak anticy-clonic waterspout and similar vortices near to, and impacting, Hong Kong International Airport, evaluating the role of horizontal shearing instability in vortex formation,   (blue/black circles) and the waterspout (red dot). Note that the Gaolan radar was not used in this study due to its distance from the storm (Map source: Google Earth).
finding similar mechanisms to those found in Wakimoto and Lew (1993), Brady and Szoke (1989), Miglietta (2019) and others. Van Den Broeke and Van Den Broeke (2015) utilised dual-polarisation observations, focusing on the predictive value of boundary layer radar observations for forecasting and documenting dual-polarisation signatures (e.g. an arc of high differential reflectivity; Kumjian and Ryzhkov, 2007) associated with the waterspouts. Based on initially weak low-level reflectivity and high differential reflectivity aloft, they inferred that stretching of pre-existing vertical vorticity along a boundary was the primary mechanism for waterspout formation.
Just after noon local time (utc + 8h) on 1 June 2021, a waterspout was sighted and photographed near the Pte. da Amizade Bridge near Ilheu Kai-Kiong (Figure 1). A report on the waterspout can be found in the online archive of the Macao Meteorological and Geophysical Bureau 1 and in the newspaper Headline News 2 , in Chinese. The novelty of this case is that the waterspout was observed by a number of phased-array, Doppler polarimetric weather radars at the same time. This weather radar network has only recently been installed, and there were three radars scanning the storm from a range of less than 40km. These Doppler weather radars enable a detailed analysis of the structure of the parent storm associated with the waterspout.
This report details the location, strength and evolution of the smallest resolvable circulation within the storm between about 0446 and 0457 utc (all times given hereafter are utc), based on an analysis of the radar data. Since the waterspout's location and path we not mapped in detail, radar observations are used to estimate its position to the degree that it can be identified in the data. Dual-Doppler and axisymmetric wind analyses and discussion of the dualpolarimetric signatures seen with the storm are presented.

Phased-array weather radar network
Three phased-array weather radars (PAWRs) installed at Hengqin, Qiao and Huangyangshan with antenna heights above sea level (ASL) of 482, 205 and 595m, respectively, and with central frequencies ranging from around 9.3 to 9.5GHz ('X-band'), are used in this analysis. In addition to beam scanning using mechanical rotation, the PAWRs scan with a multi-beam forming technique, that is, 12 electronic scans (0. 9°, 2.7°, 4.5°, 6.3°, 8.1°, 9.9°, 11.7°, 13.5°, 15.3°, 17.1°, 18.9° and 20.7°) in elevation, allowing the scanning of a complete volume every 90s. The half-power beam width in azimuth and elevation, and the range resolution are 3.6°, 1.8° and 30m, respectively. Figure 2 shows the location of the radars and the position of the circulation at 04:53:30 when it most likely produced a waterspout, based on the appearance of a small-scale couplet in the velocity data of the nearest radar (Hengqin). Four radars are shown, but the southernmost radar (Gaolan) was too distant from the storm to be included in this analysis. Of the remaining three, only one (Hengqin) is close enough, ~10km, to possibly resolve the waterspout itself. The Hengqin and Huangyangshan radar pair was used in the dual-Doppler analysis of the storm. Due to a lack of useful data near the waterspout from the Qiao radar, the Qiao and Huangyangshan radar pair was only used for verifying the vortex structure, based on the analysis from the other radar pair. The Hengqin radar was used for single Doppler analyses and analysis of dual-polarisation signatures due to its proximity to the storm.

Evolution of the storm and its circulation, as observed by the Hengqin radar
Identification of the waterspout in the velocity data is crucial in this analysis. For times prior to 04:55:00, there is only one velocity couplet that could be the waterspout.  Subsequently, however, there is a velocity couplet with two smaller couplets embedded within the larger couplet, making it a 'multiple-vortex' (mv) circulation. The waterspout could be associated with either, but the western sub-vortex is stronger, tighter and also shows a reduction of ρ hv , the crosscorrelation coefficient, indicating lofted debris (grass, paper, etc. or larger items). This study thus identifies the western subvortex to be associated with the waterspout.
With the waterspout identified in the radar data from Hengqin radar, it can be roughly tracked, as shown in Figure 3, with a horizontal positional uncertainty of about 150m. From this track, it is calculated that the waterspout had an average propagation speed of about 3ms −1 towards the northeast. Table 1 details the distance, diameter and intensity of the circulation. Intensity is taken to be the difference between the inbound and outbound maximum Doppler winds (DV) recorded in the lowest sweeps (elevation angle 0.9°) of the Hengqin radar. For times with 'mv' , derived metrics are associated with the western sub-vortex. Diameter is the distance between the maximum inbound and outbound Doppler velocities, or twice the mean radius of maximum winds (RMWs). At the range of about 10km, the smallest distance between azimuthal radar data samples is about 175m. Therefore, any waterspout smaller than about 350m (2 × 175m) is not well resolved, and DV is poorly estimated. The maximum estimated DV was 27ms −1 at 04:55:00, but this was measured across the larger vortex containing two weaker sub-vortices (shown in Figure 4). The strongest DV measured across a single sub-vortex was 21.3ms −1 at 04:53:30.
At 04:53:27, the visual appearance of the Doppler velocity signature of the subvortex, with narrow velocity extrema on opposite sides (northwest and southeast) of the vortex ( Figure 5), is consistent with the sub-vortex being tilted from the vertical direction towards either the southeast or northwest, approximately perpendicular to the radar beam from the Hengqin radar. This temporarily allows for metrics of the vortex to be resolved using the 30m-scale radial data spacing as opposed to the coarser azimuthal spacing. The data thus indicate the vortex may have been as small as ~30-60m in diameter at that time. This is much smaller than the 170m azimuthal resolution can reveal normally, when the vortex is vertical. The implication is that the strength of the vortex is underestimated due to averaging within the relatively large sampling volumes of each radar data point.

Dual Doppler analyses
The Doppler velocity data were quality controlled before further processing, removing data points that deviated by more than 2ms −1 from the median of the immediately surrounding data. The resulting gaps were filled by averaging the surrounding remaining data. Radar volumes spanned elevation angles from 0.9° to 20° about every 2°. Using a method similar to Kosiba et al. (2013), radar data from the Huangyangshan, Hengqin and Qiao radars were interpolated to a Cartesian grid using a two-pass Barnes analysis (Barnes, 1964;Majcen et al., 2008) with grid spacing and smoothing parameters (kappa) determined by the most  At this time, the azimuthal data spacing at the vortex is 160m, while the radial data spacing is just 30m. The elongated nature of the couplet suggests the vortex is slanted to some degree out of the vertical, in a direction perpendicular to the radial direction and allows a temporary opportunity to resolve the vortex in the radial direction to as fine as 30m. The sharp inbound maxima (dark green pixels) indicate that the vortex is about 30-60m in diameter.
poorly resolved radar data at the location of the circulation (the Huangyangshan radar). This step transfers data from the irregularly spaced radar grids to the regularly spaced Cartesian grid using a distancesensitive weighting function. At a range of 36km from the waterspout, the azimuthal spacing of data from the Huangyangshan radar is 630m. This resulted in the choice of Cartesian grid spacing of 100, 100 and 200m in the X, Y and Z directions, respectively, with X increasing eastward, Y increasing northward and Z increasing upward. The analysis grid origin is defined to be at the Hengqin radar.
At a given point within the gridded domain, the quality of the dual-Doppler retrieval of the full horizontal wind vector is improved as the viewing angle between the two radars increases, being best at 90°. A minimum beam-crossing angle of 30° was required. For this study, two such pairs of radars (Hengqin, Huangyangshan and Qiao, Huangyangshan) are used to get two estimates of the full wind at each point. However, only the first pair had data coverage at the location of the waterspout from both radars. Therefore, dual-Doppler analyses from the Qiao-Huangyangshan pair are presented only for comparison with the Hengqin-Huangyangshan pair. Figure 6 shows the arrangement of the Hengqin, Huangyangshan dual-Doppler lobe, the area within which the radar beam viewing angle is greater than 30°, shaded green, along with the dual-Doppler retrieval domain used in this study (red box). The waterspout is indicated by a red dot. Figure 7 depicts the evolution in time of the low-level (400m ASL) storm structure derived from Hengqin, Huangyangshan dual-DV retrievals. Vectors depict the horizontal wind. The wind vectors are dominated by the mean ambient wind blowing through the domain, almost masking the vortex. Subtraction of the mean wind at each height of the domain allows one to view the vortex more clearly. The mean wind varies with height between about 12ms −1 from the south-southeast at low levels to 12ms −1 from the west at higher levels with small changes during the analysis period. To allow easier understanding, instead of subtracting the mean wind at each height, we subtract a constant 10ms −1 northward wind for all low-level plots. This allows one to better see changes of wind with time. If desired, the full wind can be quantified by adding back a 10ms −1 northward wind at each point. Colours represent the vertical wind (W), determined by the vertical integration of the horizontal divergence, with positive being upward. It should be noted that the calculation of W depends on the presence of low-level data, which is not present in this case, especially for the more distant Huangyangshan radar, so it  should be considered a qualitatively valid estimate only. Blue lines show the vertical vorticity (spin of horizontal winds) which is derived from the dual-Doppler retrieval. Reflectivity values (black lines) are taken from the nearby Hengqin radar.

Figure 7. Dual-Doppler retrievals at Z = 400m ASL at various times. W is indicated by red colours for positive values indicating upward vertical velocities (blue for negative/downward), reflectivity (from the Hengqin radar) is indicated by black lines starting at 40dBZ, increasing in intervals of
The broad-scale storm structure is similar at all times. The reflectivity field exhibits a large 'hook-like' feature in the southwest part of the storm, which corresponds to a broad circulation, as evidenced by the horizontal winds. Another circulation is apparent in the northern part of the domain, probably associated with another storm. To the west of the circulation is an area of upward motion, which is in contrast to most dual-Doppler analyses of tornadic storms, where the 'rear-flank downdraft' usually appears (Kosiba et al., 2013). The Huangyangshan, Qiao dual-Doppler analyses show the updraft as well. Most analysis features are similar between the two dual radar pairs (see Figure 8), though one must keep in mind that both pairs employ the Hengqin radar. Figure 9 depicts the vertical variation of the winds at 04:53:30, which is representative of all the times. For clarity, the mean wind at each height (rather than a constant wind, as was done for Figure 7) has been subtracted due to its variation with height. Table 2 lists the eastward component (U) and the northward component (V) of the mean wind at each height that is displayed in Figure 9. The circulation remains discernible to at least the top of the data (4.5km ASL), beyond what is shown in the figure, becoming less well organised with height, and seems to shift some to the northeast.

Ground-Based Velocity Track Display analysis
The dual-DV retrieval was constrained by the resolution of the most distant radar. Thus, the waterspout itself could not be resolved; just the broader-scale circulation within which it was embedded. A finer spatial-scale analysis of the circulation using single-Doppler data is possible using the better spatially resolved data from the closest radar (Hengqin). Using the Ground-Based Velocity Track Display (GBVTD) algorithm (Lee et al., 1999;Lee and Wurman, 2005), which assumes an axisymmetric vortex, profiles of the tangential (V tan ) and radial (V r ) components of the wind relative to the vortex centre as a function of range from the centre and height ASL are produced. Within a radar volume (where sweeps increase in elevation covering a 3-dimensional space), a vortex centre point is chosen manually at each height. Doppler data are sampled in rings concentric with the centre, and for each ring, V tan and V r are calculated based on that data. From these two horizontal components, axisymmetric W is calculated from the vertically integrated divergence, with W = 0 as a lower boundary condition. Figure 10 shows V tan as colours, and the secondary wind circulation (V r and W) as vectors, all as functions of distance from the vortex centre, and height, at several times. The plots can be thought of as profiles of the vortex's wind, with the   Table 2) has been subtracted to more clearly show the circulations. The broad scale circulation is present at least as high as this data extends, though it becomes less well organised with height.

Table 2
Mean wind used in dual-Doppler plots of Y-axis representing the centre of the vortex as it extends upward. From 04:45:50 to 04:53:30, the results appear to be consistent in time, with the radius of maximum wind (RMW) about 500m wide near the surface, and a pronounced downdraft within the RMW at most heights. Maximum V tan gradually increases from about 3 to 7ms −1 with time. Then the vortex becomes differently organised, consisting of multiple smaller vortices (RMWs ~150m as noted in Table 1) embedded within a larger vortex during the 04:55:00 and 04:56:40 volumes. These are indicated in the Doppler radar data in Figure 4, which is annotated with a large circle representing the parent vortex, and small circles representing the sub-vortices. This type of organisation is common with many tornadoes (see Wurman and Kosiba, 2013), where small-scale circulations are embedded within larger-scale vortices and defining 'what is the tornado' is often difficult and subjective.
Since the GBVTD results are, by definition, centred on the vortex, the results are influenced by the choice of that centre, which depends on which vortex (shown in Figure 4) is being analysed. Therefore, two different GBVTD retrievals are performed for these two volumes: one centred on the larger, parent vortex and the other centred on the dominant (western) sub-vortex. The results centred on the sub-vortex show its structure -a smaller, weaker (and less-well resolved) circulation. Centring on the larger, parent vortex shows the stronger, wider circulation. Keep in mind that both of these circulations are, in reality, present simultaneously. Both results are valid but represent different vortices. Figure 11 compares the two GBVTD retrievals for these two volumes. These plots suggest that the sub-vortex extends to a height of 2km where it merges with the larger vortex. The small vortex had upward motion while the surrounding parent vortex had downward motion. It is likely the weakness of the smaller vortex is primarily due to the poor radar resolution (relative to its size) and would be considerably stronger if it was fully represented in the data. Equally, the upward motion in the smaller vortex would likely be stronger if data nearer the surface were available.
Photographs of the waterspout (Figure 1), showing a very narrow funnel (on the order of 200m in diameter) which suggest that it was associated with one of the sub-vortices rather than the much larger parent vortex.

Dual-polarisation signatures
Dual-polarisation signatures in radar data, in particular the differential reflectivity (ZDR) and ρ hv mentioned earlier, can provide additional insight into the storm kinematics (Kumjian and Ryzhkov, 2007). The type of storm can be ascertained, in part, from the presence or lack of a 'ZDR Arc' , an enhancement of ZDR north of the circulation along the inner edge of the strong reflectivity core. Its presence indicates raindrop size sorting which is indicative of shear in the environmental wind field supportive of the formation of supercells, which are generally capable of producing stronger tornadoes than non-supercell storms. This kind of information can be useful in making the decision to issue a tornado warning. Of possibly greater interest is that the presence of a tornadic circulation on the ground can be inferred by a co-located reduction of ρ hv due to the lofting of even light debris, which tumbles in a random manner, reducing ρ hv . Figure 12 shows these two dual-polarisation variables, ZDR and ρ hv , along with Doppler velocity (VEL) and reflectivity (Z) at the same time, as observed by the Hengqin radar. A white circle marks the location of the waterspout in each.
The ZDR arc in this case generally mirrors the shape of the high reflectivity contour rather than being concentrated along the edge of the strong reflectivity just north of the waterspout and does not evolve much during the analysis period. Thus, it may not be indicating much about the storm, or it may be indicating weak shear in the environmental wind, as evidenced by the small variation with height of the mean wind noted in the dual-Doppler retrievals ( Table 2). The persistent hook shape in the reflectivity field with the persistent co-located deep circulation indicated in the dual-Doppler-derived wind field is a more compelling indication that this storm is likely to have super-cellular characteristics.
The ρ hv field shows a rather pronounced drop at the location of the waterspout from 04:53:30 until the end of the analysis period. This is an indication that the waterspout lofted some (probably light) debris. The lofting of water droplets alone is not sufficient to reduce ρ hv since water droplets do not exhibit tumbling-related radar return variations (Kosiba et al., 2012). Whether or not this implies that the circulation touched land is unknown, as inflow may have carried light debris into the waterspout from just offshore, or from the bridge.
This event helps to verify that these radars are capable of detecting even weak tornadoes using their various products, including dual-polarisation capabilities.

Conclusions
A waterspout-producing storm -likely a supercell -at Macao, China, is analysed using dual-and single-Doppler analyses. Dual-Doppler retrievals show a persistent, large-scale circulation within the storm extending through the depth of the data (4.5km ASL), but do not resolve the smaller waterspout due to its distance from the radars. The Doppler velocity signature of the waterspout is evident to some degree in data from the closest radar only. Initially broad and diffuse, the circulation began to show sub-vortices, though initially it is dif-ficult to discern which vortex was associated with the observed waterspout. A ρ hv drop evident in the stronger, tighter sub-vortex identified it as likely being the one associated with the waterspout. It also indicated that the waterspout had ingested light debris from on shore, but of what and how much is not known.
The vortex structures are analysed further using the GBVTD single-Doppler axisymmetric retrieval method. GBVTD analyses centred on the large and on the sub-vortex indicated that the sub-vortex extended to a height of 2km before becoming indistinguishable from the larger vortex within which it was embedded. The small vortex, whose magnitude was underestimated due to poor resolution, had upward motion (probably also underestimated) while the larger one has downward motion.
Dual-polarisation data from the Hengqin radar, along with the deep persistent circulation evident in the dual-Doppler retrievals suggest that this storm was a weak supercell.
The data collected during this event show the usefulness of southern China's phasedarray radar network in detecting waterspouts in an operational setting, and for providing data for research, although the usefulness of the research was somewhat limited by the height of the radars above the sea-surface, rugged terrain and distance between the radars on this occasion. When events like this one occur in the future, it is hoped that similar operational data will be analysed to confirm and extend these results. It is possible that more favourable geometry will permit observations at lower levels.