## Abstract

Formulas are obtained for observed circulation around and contraction rate of a Doppler radar grid cell within a surface of constant launch angle. The cell values near unresolved axisymmetric vortices vary greatly with beam-to-flow angle. To obtain reliable standard measures of vortex strength we bilinearly interpolate data to points on circles of specified radii concentric with circulation centers and compute the Doppler circulations around and the areal contraction rates of these circles from the field of mean Doppler velocities. These parameters are proposed for detection of strong tornadoes and mesocyclonic winds. The circulation and mean convergence around the Union City, Oklahoma, tornado of 24 May 1973 are computed. After doubling to compensate for the unobserved wind component, the circulation (1.1 × 10^{5} m^{2} s^{−1}) agrees with a previous photogrammetric measurement. The mature tornado was embedded in a region, 6 km in diameter, of nearly uniform strong convergence (~5.5 × 10^{−3} s^{−1}) without a simultaneous mesocyclone. A model of a convergent vortex inputted to a Doppler radar emulator reproduces these results. Moving the model vortex shows that for a WSR-88D with superresolution, the circulation is relatively insensitive to range and azimuth. WSR-88D data of the 31 May 2013 El Reno storm are also analyzed. The tornado formed in a two-celled mesocyclone with strong inflow 5 km away. In the next 8 min the circulation near the axis doubled and the areal contraction rate at 5 km increased by 50%. This signified a large probability of strong tornadoes embedded in powerful storm-scale winds.

## 1. Introduction

Researchers first used pulsed Doppler weather radars in the 1970s. They quickly realized that characteristic patterns in the Doppler velocity field observed with a single Doppler radar revealed the presence of mesocyclones and tornadoes (Brown et al. 1978; Burgess et al. 1993). Attention turned in the 1990s to automatic detection of mesocyclones and tornadoes. Meteorologists and automated algorithms customarily measure the strength of a convergent vortex by the velocity difference (delta-*V* or Δ*V*) between the two peaks in a characteristic velocity couplet (e.g., Mitchell et al. 1998; Stumpf et al. 1998; Kuster et al. 2015). Closely related parameters are the rotational velocity *V*_{rot} = Δ*V*/2 and the shear *S*_{D} = Δ*V*/*D* where *D* is the distance between the peaks. For an unresolved vortex of given strength these measures vary considerably with range from the radar. More complicated parameters are the observed or Doppler circulation Γ_{ob} around a circle of specified radius (Davies-Jones and Stumpf 1997), the Doppler circulation of the velocity couplet Γ_{D} (Davies-Jones and Wood 2006), and the excess rotational kinetic energy (ERKE; Donaldson and Desrochers 1990). At a single height, the rotational kinetic energy (RKE) is proportional to $\Gamma D2$, and so is similar to, but more changeable than Γ_{D}. The ERKE is the RKE calculated by replacing Γ_{D} with Γ_{D} − Γ_{M} where Γ_{M} is a “threshold circulation” associated with a minimal mesocyclone of the same apparent radius as the observed one. If Γ_{D} − Γ_{M} is negative, the ERKE is set to zero, which is thus the minimum possible value. To test the utility of ERKE compared to Γ_{D}, Davies-Jones and Wood (2006) used a virtual Doppler radar to “observe” an analytical contracting vortex flow (Rott 1958). They found that ERKE remained zero during the initial spin up of the vortex and increased very rapidly as the mesocyclone contracted to a tornado. Thus, thresholding and volatility acts to decrease the lead time of a hypothetical warning based on ERKE compared to one based on Γ_{D}.

Another approach involves estimating azimuthal shear by fitting a least squares plane to the Doppler velocity field in a local grid of adjacent discrete azimuths and range bins (Mahalik et al. 2019). This method yields the neighborhood’s mean azimuthal shear and a mean along-azimuth shear. The data are first subjected to a median filter, which for each 3 × 3 subset of grid points finds the median value and assigns it to the central point. This method was tested with a collection of analytical vortices of various strengths and radii from 1 to 8 km as observed with a WSR-88D simulator. Because observed azimuthal shear is range dependent, estimates of the rotational intensity was good only for ranges up to 90 km (probably less for the smallest vortices in the group). The method is unsuited for estimating the strengths of unresolved vortices because of data smoothing and use of derivatives of Doppler velocity.

In contrast to Δ*V* and azimuthal shear, Davies-Jones and Wood (2006) found that the circulation of a Doppler velocity signature Γ_{D} provides a fairly good estimate of the actual value for the vortex, and is relatively insensitive to range, beamwidth, and stage of vortex evolution. Rapid contraction of a rotation signature with large circulation provides indication of imminent tornadogenesis to a weather forecaster (Kuster et al. 2015).

Circulation Γ_{ob} is an important parameter for detection of potential and actual tornadoes because it is a far-field measure of potential vortex strength and unlike the other parameters discussed above is independent of where the apparent velocity peaks of an unresolved vortex happen to lie. The far field of a tornado naturally approximates a potential vortex with constant circulation because the precise flow structure of the core flow becomes inconsequential rapidly with distance. Based on a proposal by E. N. Rasmussen (1995, personal communication), Davies-Jones and Stumpf (1997) tested circulation around a circle 2 km in radius as a tornado-detection parameter. Consider, for example, the circulation around a circle centered on the axis of a Rankine combined vortex (RCV), which crudely represents a tornado or mesocyclone. The circulation Γ in the far-field (outer potential-vortex) part of an RCV is constant and is equal to 2*π* times the radius *ρ* of the circle times the tangential velocity *υ* around the circle. Hence a measured circulation of 6.28 × 10^{4} m^{2} s^{−1} around a circle of 1500 m might be indicative of an intense tornado with a maximum tangential wind of 100 m s^{−1} at a radial distance of 100 m. Alternatively, it might be symptomatic of (i) a much wider and weaker vortex, which could spin up into a powerful tornado, or (ii) multiple vortices within the circle. A much smaller circulation around the circle, say 6.28 × 10^{3} m^{2} s^{−1}, is enough for only a narrow tornado. The circulation of a mesocyclone is typically 4 × 10^{5} m^{2} s^{−1}, based on a maximum tangential velocity of 21 m s^{−1} at a radial distance of 3 km. The contraction rate of the circle is also an important parameter in the spinning up of vortices.

We can compute circulation around and contraction rate of several specified circles quickly and easily from operational Doppler radar data. The measurable radar coordinates (*r*, *α*, *β*) are slant range *r* (arclength distance along a stationary ray from the radar to a measurement point), ray azimuth angle *β* measured clockwise from due north, and the launch angle *α* (elevation angle at the radar) of a ray. During an operational volume scan of a WSR-88D the velocity data are recorded on a 3D grid with the data collected sequentially along radials at constant range-gate spacing Δ*r*, then at successive azimuth angles, and finally at different elevation angles specified by the particular volume coverage pattern (VCP) in use. Therefore, the fastest way to detect the presence of dangerous vortices is to determine kinematic quantities from the data in surfaces of revolution of constant launch angle (hereafter *α* surfaces). Since a single Doppler radar only observes flow along the curved beam, we can compute relevant kinematic quantities such as circulation only partially. We will refer to the Doppler-observed parts of circulation as Doppler circulation, and similarly for areal contraction rate, etc.

To compute the height of data points, we account for the bending of rays by atmospheric refraction by assuming a standard atmospheric stratification and ray curvature that is constant along rays but varies with the cosine of the launch angle (Davies-Jones et al. 2019, hereafter DJWA). This is done in DJWA for different earths, namely, the actual Earth, an enlarged equivalent earth on which the rays are straight, and a flat earth with modified ray curvature. On all these earths, height above radar level is the same function of slant range (to an accuracy of a few meters at most) provided that the planetary curvature minus the curvature of a ray launched horizontally is held constant (Fig. 1).

Kinematic quantities within an *α* surface such as shear, rotational velocity, convergence, and circulation, depend on the ray curvature but not on the planetary curvature. Over the domain of weather radar, the ray bending is negligible under standard refraction so we may neglect ray curvature for the Earth. Thus, we may assume that (*r*, *α*, *β*) are spherical coordinates when computing the above kinematic quantities.

Ideally a measure of vortex intensity should be fairly independent of the range and azimuth of the vortex from the radar. Since the intensity of a long-track tornado varies along its path, we determine the range dependence of various vortex-strength measures using a simple analytical vortex flow and a Doppler radar simulator (as in Davies-Jones and Wood 2006). By moving the axis of the vortex, we can generate the field of mean Doppler velocities at different ranges. The analytical vortex flow is over a flat earth so we employ the DJWA flat-earth model, which compensates for the lack of earth curvature with a precisely modified ray curvature. With this curvature we obtain the flow as a function of radar coordinates as follows. For each radar location we locate the points (*X*, *Y*, *H*) in the flat-earth Cartesian geometry corresponding to each radar grid point (*r*_{i}, *α*, *β*_{k}) on a particular *α* surface. We compute the flat-earth ray slope angle and the Doppler component of velocity at these points. We then use the radar emulator described in Davies-Jones and Wood (2006, 1046–1047) to compute the mean Doppler velocity at the radar grid points.

The present paper presents a methodology for computing the circulation and areal contraction rates of vortex flows from radar observations. The proposed methodology is applied to both actual and simulated observations. Section 2 derives the necessary formulas. As one example, we determine observed circulation and contraction rate of the 1973 Union City, Oklahoma, tornado (Brown et al. 1978) in section 3 and use a simple vortex model together with the above radar emulator to reproduce these observed quantities in section 4. In section 5 we move the model vortex to investigate how the observed circulation varies with range and azimuth. Section 6 describes a more recent example, the 31 May 2013 El Reno, Oklahoma, tornado, observed with a WSR-88D. The final section restates our main points.

## 2. Mathematical formulation

We assume that Earth’s surface is a sphere with radius *a*_{⊕} (≈6371 km) equal to the distance from the center of Earth to the radar and that the ray curvature *κ* owing to atmospheric refraction is constant along rays but varies from ray to ray according to the cosine of the launch angle. We use Cartesian coordinates (*X*, *Y*, *H*) associated with the tangent plane at the radar as well as the radar coordinates (*r*, *α*, *β*). The origins of both systems are at the radar antenna. Angles are measured in radians unless stated otherwise. The position vector in the Cartesian system is

where **i**, **j**, and **k** form a right-handed orthonormal basis with **i** eastward, **j** northward, and **k** upward at the radar antenna. These constant vectors are either parallel or normal to the tangent plane at the radar and do not follow Earth’s surface.

The transformation **X** = **T**(*r*, *α*, *β*) to Cartesian coordinates from radar coordinates is

where

by trigonometric identities and

is the ray curvature and *κ*_{0} is the curvature of rays launched at *α* = 0° [see Eqs. (2)–(4) of DJWA]. For standard refraction on the actual Earth *κ*_{0} = (5.76 *a*_{⊕})^{−1}. The corresponding modified curvature for use on a flat earth is −(1.21 *a*_{⊕})^{−1} (Fig. 1). From (3) and (4) it follows that

which verifies that the stationary rays (curves of constant *α*) are circles with centers at Σ = tan*α*/*κ*_{0}, *H* = −1/*κ*_{0} and radii sec*α*/*κ*_{0}. Note also that

At this point we can proceed without further approximation by regarding (*r*, *α*, *β*) as nonorthogonal coordinates (Margenau and Murphy 1956, 192–197) and retaining the appropriate ray curvature for the earth. The computer program used to generate the results in this paper in fact does just this. Alternatively, with negligible loss of accuracy as long as the planetary curvature is set to the equivalent-earth one (Fig. 1), we may set the ray curvature to zero, which is equivalent to assuming that (*r*, *α*, *β*) are spherical coordinates. We use this assumption here to simplify the following mathematical description. The formulas derived using nonorthogonal coordinates reduce to those derived below in the limit of vanishing ray curvature. To justify omission of the ray curvature, we expand the trigonometric function in (6) into series in the small quantity *κ*_{0}*r* (which for a long range of 225 km is equal to 0.006 for standard refraction on the actual Earth and −0.03 for the equivalent refraction on a flat earth). This yields

The difference between the curvilinear and straight-line slant ranges is less than 0.4 (9) meters for *r* = 225 km in the real (flat) Earth case. Therefore, for computing the kinematic quantities in a surface of constant launch angle we can neglect curvature over the ranges relevant to weather radar. In contrast, an observation point’s height above radar level is sensitive to the difference between ray curvature and planetary curvature. Height and ground range relative to the planetary surface and slope angle of the ray relative to the local horizontal are computed here as in DJWA.

The circulation Γ of the radar targets around a circuit *C* that encloses an area *A* in an *α* surface (*dα* = 0) is by definition

and the observed circulation is

where *V*_{D} is the mean Doppler velocity. The mean vorticity is the circulation divided by the area.

Circulation is calculated from the gridded mean Doppler velocity data. The grid points in a surface of constant launch angle are located at *r* = *r*_{i}, *β* = *β*_{k}, where subscripts *i* and *k* are used to label the discrete ranges and azimuths and *r*_{i} = *i*Δ*r*. The *ik* grid cell is defined as the area in the surface enclosed by successive slant-range circles *r*_{i} and *r*_{i+1} and successive radials *β*_{k} and *β*_{k+1}.

We evaluate a general line integral

in the *α* surface as follows. Here *f* and *g* are general scalar variables. We define a counterclockwise sequence of points *P*_{m} = [*r*_{m}, *α*, *β*_{m}] with cyclic labeling *m* = 1, 2, …, 2*M*, around *C*. These points are either grid points or interpolation points. If *P*_{m} is an interpolation point, we use bilinear interpolation to obtain values there from the gridded data. Assuming that *f* varies linearly with *g* between adjacent points, the contribution to the line integral from the curve segment *P*_{m−1}*P*_{m} is

Summing all the contributions yields

where

(Davies-Jones 1993). Applying (12) to (9) gives us for the observed circulation

In particular, the observed circulation around the *ik* grid cell is

Meteorologists can detect mesocyclones from a plan position indicator (PPI) display with the cells outlined and the circulation values of each cell displayed. We can easily calculate the observed circulation around a group of cells because circulations (and line integrals in general) are additive. The circulation around the boundary of any area that is partitioned into subareas is equal to the sum of the circulations around the perimeter of each subarea because the line integrals along “interior sides” (ones that are shared by two of the subareas) cancel (Fig. 2; Petterssen 1956, p. 127). Note also that the observed grid cell circulations near an unresolved vortex suffer from a confusing spurious quadrupole pattern as described below.

To find a formula for areal contraction rate on an *α* surface we start with a line integral for the area *A*, which we find as follows. On an *α*-coordinate surface, the vector surface element is

where **u**_{α} is the unit vector in the direction of increasing *α*. The component in the **u**_{α} direction of the curl of a continuous and differentiable vector field **B** is

in spherical coordinates where *B*_{r} and *B*_{β} are the radial and azimuthal components of **B**, respectively. By choosing *B*_{r} = *βr* cos*α* and *B*_{β} = 0, we obtain an integral equal to *A*, namely,

We then use Stokes’s theorem to obtain the line-integral formula for area

where *A* is positive when *C* is traversed counterclockwise. Via (12) the computational formula is

The areal contraction rate^{1} is the rate −*dA*/*dt* at which the area *A* is contracting materially within the *α* surface and the mean convergence is the areal contraction rate divided by the area. They are important because contraction stretches and amplifies the vorticity component normal to this surface. Contraction of a circle circumscribing a mesocyclone in an *α* surface would indicate that the mesocyclone is spinning up. Note that the contacting area is semimaterialistic because it does not move with the flow normal to the *α* surface. Thus, in the following analysis we must regard *α* as constant and set the surface-normal velocity component to zero.

Applying the −*d*/*dt* operator to (19) with the above proviso yields

where *r* cos*α*(*dβ*/*dt*) is the azimuthal velocity component *V*_{β}. Now

With this result (21) becomes

after integrating by parts. The terms on the right side are, respectively, the observed and unobserved fluxes along the *α* surface into *A*. By the line-integral method, the observed areal contraction rate (with *V*_{β} set to zero) is

where the function Λ is defined in (13). Like circulations, contraction rates are additive. The areal contraction rate of a union of contiguous subareas is the sum of the areal contraction rates of each subarea.

Except on rare occasions, a tornado is too small to be resolved by a Doppler radar. The radar usually observes just the circulating winds outside the vortex core where the flow resembles a potential vortex. The Doppler vorticity and divergence fields in this region contain spurious “quadrupole” patterns (Fig. 3 and appendix) because the radar measures just the along-beam velocity component of this flow. The spurious quadrupole contributions to the observed circulation around and contraction rate of a circle centered on the vortex average out. Thus, the best choice for *C* on a given *α* surface is a circle. Furthermore, computing the observed circulation around a circle of radius *ρ* in the constant-*α* surface provides a more precise standardized measure of circulation for indicating the trend of a moving unresolved vortex’s intensity or for comparing the strengths of different unresolved vortices because it is fairly insensitive to range and azimuth. For an accurate assessment of circulation, *ρ* should be greater than the azimuthal spacing (in units of length) at the range of the vortex.

To construct the circle *C* of radius *ρ*, we first estimate the central point (*r*_{0}, *α*, *β*_{0}) where the axis of the vortex (mesocyclone or tornado) intersects the surface of constant *α* and let (*r*_{c}, *α*, *β*_{c}) represent points on *C*. In the spherical coordinates the circle is the intersection of the sphere of center (*r*_{0}, *α*, *β*_{0}) and radius *ρ* with the surface *α* = constant. The equation of this circle is

This simplifies to

after use of a trigonometric identity. Thus, the azimuths of the circle points as a function of their slant ranges are

The circle points that lie on the stationary ray through its center are

The circle points that lie at the same slant range as the center are

To a good approximation these points are the azimuths of the rays that are tangents to the circle. We use (28) and (29) to determine whether a given circle fits into the analysis region.

To compute the observed circulation around and the areal contraction rate of a circle *C* centered on a mesocyclone or tornadic vortex signature (TVS) in an *α* surface we used the following procedure. The sequence

defines uniformly spaced slant ranges spanning the spread of slant ranges across the circle [see (28) and Fig. 3]. We use (27) to compute the azimuths $\beta m+$ and $\beta m\u2212$ at which the slant ranges intersect or touch the circle (Fig. 3). At each of the points (*r*_{m}, $\beta m+$) and (*r*_{m}, $\beta m\u2212$) on *C* (2*M* points in total) we obtain the Doppler velocity *V*_{D} by locating the grid cell that the point lies in and using bilinear interpolation. We then order the points in a counterclockwise sequence around the circle and use (14) and (24) to compute the observed circulation and areal contraction rate. We can vary the radius *ρ* of the circle through separate implementations of the algorithm. Because circulation is additive, the circulation around a larger circle is the sum of that around a smaller concentric circle plus the circulation around the annulus between the circles (as explained in Fig. 2). Areal contraction rate also has this additive property. The mean convergence within a circle is its areal contraction rate divided by *πρ*^{2}, its area.

## 3. Circulation and areal contraction rate of the 1973 Union City tornado

To illustrate our procedure, we analyzed data obtained by the Norman Doppler radar at 1546 CST in the tornadic storm that struck Union City on 24 May 1973 (Lemon et al. 1978; Brown et al. 1978). From photogrammetry by Golden and Purcell (1978a) the maximum measured circulation of the violent tornado was 1.05 × 10^{5} m^{2} s^{−1} at 200 m radius and 90 m AGL. Radar detected the tornado as a TVS. We use the field of dealiased mean Doppler velocities (relative to the TVS motion of 10 m s^{−1} from 283°) in the 3.8° launch-angle surface (Fig. 4; Brown et al. 1978) obtained during the tornado’s mature stage (Golden and Purcell 1978b) to compute the circulation around and mean divergence within various circles centered between velocity peaks (at 51 km, 292.5°, 3.5 km AGL) in this *α* surface. At this time and height, the peak minimum and maximum velocities were *V*_{min} = −21 m s^{−1} and *V*_{max} = 36 m s^{−1} with a rotational velocity *V*_{rot} = (*V*_{max} − *V*_{min})/2 = 28.5 m s^{−1} and a peak ratio |*V*_{min}|/*V*_{max} =0.583. Prior to this observation time the tornado had damaged only a few structures (Davies-Jones et al. 1978) so probably the mean Doppler velocities at 3.5 km AGL arise mostly from the air motion with only minor contributions from the centrifuged motion of large debris, which can cause large errors in measured convergence. The recorded data have an azimuthal spacing Δ*β* of 1° and the distance Δ*r* in the radial direction between grid points is 600 m (four pulse depths). The data missing at two points were filled in using a variational technique that minimizes the gradient of Doppler velocity without modifying data at neighboring points. The 0°-launch-angle data were not used because at very low elevations the TVS was on the edge of the radar echo.

The field of cell circulations (Fig. 5) is characteristic of an unresolved RCV. The circulation field has a peak at the TVS location as expected for a concentrated cyclonic vortex core. Cells wholly in an irrotational outer flow of an RCV should have no circulation. Because a single Doppler radar does not measure the azimuthal velocity component, the peak is surrounded by an artificial pattern of cell Doppler circulations with false positive Doppler vorticity (cell circulation divided by cell area^{2}) in the far and near quadrants and false negative Doppler vorticity in the left and right quadrants (Fig. 3 and appendix). The potential vortex also causes an artificial skewed quadrupole pattern in the Doppler convergence (cell contraction rate divided by cell area) field with false Doppler divergence in the far left and near right quadrants and false Doppler convergence in the far right and near left quadrants (Fig. 3 and appendix). The false Doppler vorticity and divergence satisfy an inverse square law (i.e., decrease with the square of distance from the vortex axis).

Figure 6 shows the computed circulations around circles of various radii centered on the TVS. Results for circles smaller than 1 km in radius were unreliable because of the limited resolution of the data. The largest circle tried in this case was 3 km in radius. The algorithm automatically reduced its radius to 2.88 km to fit in the domain of Fig. 4. The circulation is practically constant for all the circles. Thus, at this time in the tornado’s life there is no evidence for the existence of a mesocyclonic vortex [for typical evolution see first figure in Lemon and Doswell (1979)]. After doubling to compensate for the unobserved wind components being set to zero ( appendix), the circulation values at 3.5 km AGL (~1.1 × 10^{5} m^{2} s^{−1}) compare favorably with the photogrammetric measurement within a few hundred meters of the ground. Note that the circulation of a vortex is constant along its length by Kelvin’s circulation theorem so measurements at different heights should be roughly the same.

The mean convergence, −(1/*A*)*dA*/*dt*, also does not vary appreciably as the radius varies from 1 to 3 km, indicating that the mature tornado was embedded in a region of fairly uniform strong convergence (~5.5 × 10^{−3} s^{−1} after doubling) that was around 6 km in diameter (Fig. 7). The tornado appears to be a contraction and spinup of a parent mesocyclone that existed prior to this time (Lemon et al. 1978). For a circle concentric with the axis of an axisymmetric flow, the ratio of circulation to areal contraction rate is a measure of tan *χ* = −*υ*/*u* where *χ* is the inflow angle, and *u* and *υ* are the average radial and tangential wind components, respectively, on the circle. The calculated inflow angle at various radial distances around the Union City TVS is shown in Fig. 8. The inflow angle increases progressively with decreasing distance from roughly 40° at 2.5 km to 80° at 1 km from the TVS. This is again a signature of a strong concentrated vortex forming in a broad updraft. Based on insights gleamed from visual and WSR-57 radar observations, Ward (1972) built a vortex chamber that successively generated tornado-like vortices within a *wide* updraft. Most chambers at the time featured a small exhaust hole and narrow updraft with nearly constant inflow angle away from the sink. As pointed out by Ward and by Davies-Jones (1976) and confirmed by the present results, there is no evidence for a narrow sink flow aloft in a tornado. (The only region in a tornado where the flow resembles a sink is in a small corner region close to the ground where the boundary layer flow erupts upward in a high-speed axial jet.)

## 4. Simulated TVS data

We simulated the Doppler velocity field that a Doppler radar would see as it scanned an idealized vortex in a convergent axisymmetric flow. The simulator (Davies-Jones and Wood 2006) assumes uniform reflectivity within a sampling volume, an effective Earth’s radius of 1.21 times the actual value, and Gaussian weighting functions in range, azimuth (with an effective beamwidth), and elevation (Doviak and Zrnić 1993). We have modified the simulator using formulas derived herein and in DJWA so that we can more accurately emulate observations by a Doppler radar on a flat earth (for which we cannot use straight rays). The virtual radar was given a half-power beamwidth of 0.81° and a range depth of 150 m to match the Norman Doppler radar. From the first figure of Brown et al. (2002) we deduce that the effective half-power beamwidth is 1.35° for the above beamwidth and a 1° azimuthal sampling interval as used in the Union City data collection. We applied the simulator to the Union City tornado using a convergent RCV flow with a maximum tangential wind 80 m s^{−1} at a radius of 220 m (Fig. 9). The corresponding circulation is 1.05 × 10^{5} m^{2} s^{−1}. The TVS was located at 292.25°, 51 km to yield a peak ratio (0.572) that is close to the observed one. The rotational velocity, 23.2 m s^{−1}, is underestimated by 18%, indicating slightly too much smoothing of the true tangential velocity profile. The axisymmetric inflow had a similar profile to the tangential wind, with maximum inflow of 8.25 m s^{−1} at 3 km and divergence of −5.5 × 10^{−3} s^{−1} inside 3 km and 0 outside (Fig. 9). The vertical wind was found by integrating the continuity equation, using a density-scale height of 10 km.

Figures 6–8 show the simulated circulations, mean convergences, and inflow angles for different sized circles in the 3.8° constant-launch-angle surface. Once again, the values of the computed circulations and mean convergences are one-half the true values because of the zero contributions from the unobserved part of the wind. The similarities of the simulated curves to the observed ones in Figs. 6–8 indicate that the simulated TVS replicates the Union City TVS well.

## 5. Effect of range on simulated TVS

To find out how measures of the idealized Union City tornado (Fig. 9) would change with range and azimuth on an operational WSR-88D, we varied the range from 25 to 225 km in increments of 25 km and set the launch angle to the lowest one in routine use, 0.5°. Typical parameter values for a WSR-88D are as follows. The beamwidth is 0.89°, the range width is 274.5 m and the range-gate spacing is 250 m. Under standard refraction the centerline height of a beam launched at 0.5° increases from 260 m to 5.2 km in the above range interval (Fig. 10). We do not consider ranges greater than 225 km because the beam becomes too wide and high. For launch angles *α* < 1.8° (≥1.8°) super (ordinary) resolution is used currently. Since we are using *α* = 0.5° in this set of experiments we assumed superresolution. Super (ordinary) resolution has an azimuthal sampling interval (ASI) of 0.5° (1°) and effective beamwidth of 1.02° (1.39°); see Brown et al. 2002. Similar to Wood and Brown (1997), we ran two series of experiments, one with the vortex axis midway in an ASI and the other with the axis at the interface between two sampling intervals. In the former (latter) the peaks of a TVS are separated by one (two) azimuthal sampling intervals.

Figure 11 shows the simulated rotational velocity for the RCV of Fig. 9. Because it varies with the linear distance between the velocity peaks in the grid, rotational velocity decreases greatly with range (Fig. 11). The data points marked by blue and red dots are for case M where the vortex axis is at the midpoint of an ASI and for case E where the axis is at an endpoint, respectively. At moderate ranges the rotational velocity varies by around 20% owing to location of the maximum tangential winds within the ASI. In case M(E) the negative peak and the positive peak in Doppler velocity are separated by an odd (even) number of ASIs, respectively. At all the case E data points the peaks are separated by two intervals and the radar signatures are TVSs. Apart from the data points at 25 and 50 km, the peaks at case M data points are separated by one interval and so the signatures are again TVSs. At the two exceptional points the peaks are separated by three intervals and the signatures are thus tornado signatures (Brown et al. 2002) for different reasons. At 25 km range the angular core diameter of the actual vortex is 1° or twice an ASI. Thus, in case E, the peaks in Doppler velocity are found one grid point either side of the grid point where the axis is located. In case M where the axis is midway between grid points, the peaks cannot be one interval from the axis because there are no grid points there. Owing to imperfect azimuthal resolution the apparent core diameter is larger than the true one (Fig. 9.26 in Doviak and Zrnić 1993). Thus, the peaks move farther away from the axis to the next grid point and the signature becomes a tornado signature (TS) as the peaks are more than two ASI apart. At 50 km range the angular core diameter of the vortex is 0.5° so the peaks would lie at contiguous grid points in case M if the beam were infinitely narrow. Since the true angular core diameter is now one-half the EBW (1.02°), an effect, the inclusion in the average of negative Doppler velocities at points to the left of the axis, acts to reduce the mean Doppler velocity at the grid point (call it 1) immediately to the right of the axis. At the next grid point farther to the right (say point 2) a second effect, the fact that the peak velocity has a smaller weight in the weighted average, acts to make the Doppler velocity smaller at 2 than at 1. However, the first effect has the opposite impact because it is less detrimental at 2 than at 1. Thus, a TS can arise when the first effect is larger than the second effect, as is the case with the TS at 50 km range.

The simulated signature is slightly asymmetric owing to the flow’s radial inflow. The peaks are at the same range but the ratio of the negative peak to positive peak velocity varies from 1 at very short range to around 0.78 at 225 km with most of this decline occurring at ranges greater than 150 km (Fig. 12). This asymmetry vanishes when the inflow is set to zero. The existence of signature asymmetry associated with axisymmetric inflow signifies that it is difficult to estimate the location of the axis within a sampling interval based just on the ratio of peak velocities. The rotational velocity seems to be unaffected by the axisymmetric inflow when the velocity peaks are at the same range.

We now demonstrate that Doppler circulation around a circle of suitable radius is a more reliable measure of vortex strength than rotational velocity because, over a large span of ranges, it is relatively insensitive to range and to the location of vortex axis within an azimuthal sampling interval (Fig. 13) since circulation is constant outside the core of the RCV. For example, with superresolution the Doppler circulation around a circle of 2.5 km radius is nearly constant for ranges up to 225 km and is almost equal to one-half the actual circulation of the idealized vortex (1.106 × 10^{5} m^{2} s^{−1}). The range span over which the circulation around smaller circles is almost constant is smaller in proportion to the radius.

To estimate the range where azimuthal resolution begins to deteriorate Doppler circulation as a measure of vortex strength we plotted twice the Doppler circulation divided by the actual circulation (the normalized Doppler circulation) as a function of *r*/*ρ* for circles of various radii *ρ* in the case where the vortex axis is at the midpoint M of an ASI (Fig. 14). The normalized Doppler circulation is close to its perfect value of 1 for very low values of *r*/*ρ* and decreases below 0.8 as *r*/*ρ* surpasses 90. Thus, twice the Doppler circulation around a circle of radius 2.5 km departs less than 20% from its true value over ranges up to 225 km In the case where the vortex axis is at an endpoint E of an ASI, the departure is only 10% (Fig. 15). For *r*/*ρ* ≤ 90 the least accurate circulation occurs when the axis is at the midpoint of an ASI. Since the circulation around a 2.5 km radius circle may be due to either a mesocyclone or a tornado, we recommend also computing Doppler circulations around a smaller circle of radius 1.5 km for detecting possible tornadoes. Twice the Doppler circulation around this circle departs by 20% from the true circulation at a range of 135 km.

Compared to rotational velocity, circulation varies less with position of the vortex axis within an ASI. In the recommended range *r* < 90*ρ* the variation in the measured Doppler circulation owing to azimuthal location is at most 11% (cf. Figs. 14 and 15).

With superresolution the Doppler areal contraction rates of circles of various radii *ρ* in the 0.5° launch-angle surface are nearly constant for ranges up to 225 km (Fig. 16) regardless of axis location. The Doppler contraction rates are practically equal to one-half times the actual convergence (5.5 × 10^{−3} s^{−1}) times the circle area *A*. The areal contraction rate of the convergent vortex flow is estimated far better than the circulation for two reasons. First, unlike vorticity, the convergence of the model flow is constant over a large radial distance (3 km) from the axis. Second, at all ranges greater than 30 km the azimuthal grid spacing, which varies with range, exceeds the constant range-gate spacing of 250 m. Thus, the Doppler convergence generally is resolved far more finely than the Doppler vorticity.

## 6. Circulation and areal contraction rate of the 2013 El Reno tornado

As a second example we analyze some dealiased data from the well observed and much studied long-lived tornado that occurred near El Reno on 31 May 2013 (Wurman et al. 2014; Bluestein et al. 2015). Wakimoto et al. (2016) present a detailed map of the tornado damage path. We used the data from the KOUN WSR-88D at two different times, 2303 UTC, which is a minute or less after the beginning of the tornado (Kuster et al. 2015; Seimon et al. 2016; Bluestein et al. 2019), and 2311 UTC when the tornado was well established. A TVS with height continuity did not appear in the KOUN data until around 2305 UTC. The range-gate spacing is 250 m and the azimuth increment is 0.5°. The launch angles *α* at 2303 and 2311 UTC were 0.52° and 0.97°, respectively. Figure 17 shows the radar reflectivity and Doppler velocity fields at these times and elevation angles and circles *ρ* = 1, 2, …, 5 km centered on the signatures in the constant-*α* surfaces. Relative to KOUN the center of the mesocyclone signature at 2303 UTC was situated at 296.4° azimuth, range 62 km and the central beam height was 812 m above radar level. The corresponding values for the TVS at 2311 UTC were 295.75°, 55.5 km, and 1170 m, respectively. The El Reno storm is at a slightly longer range that the Union City TVS.

The steady increase of Doppler circulation from *ρ* = 1 to 3.5 km (Fig. 18) at 2303 UTC indicates a mesocyclone with an observed circulation value of 1.7 × 10^{5} m^{2} s^{−1} or 3 times the observed circulation of the Union City tornado. At this time the negative (positive) areal contraction rate of circles with radii less (greater) than 2.5 km (Fig. 19) implies that the mesocyclone was a two-celled vortex with a central downdraft (Trapp 2000). Therefore, the tornado probably formed away from the mesocyclone center (as in the first figure of Lemon and Doswell 1979). The zone of positive areal contraction rate extends from *ρ* = 2.5 km to past *ρ* = 5 km. The areal contraction rate at *ρ* = 3 km is 9 times the corresponding Union City value. However, the Union City data were at a higher elevation, 3.5 km, which makes comparison between the cases of convergence-related data problematic. Clearly the outer cell of the El Reno mesocyclone at 2303 UTC is embedded in strongly convergent flow. We define the average tangential velocity around each circle, *V*_{tang}(*ρ*), as 2Γ_{ob}/2*πρ*, where the leading factor of 2 factors in the contribution from the unobserved wind. This variable has a flat maximum of 14.8 m s^{−1} near *ρ* = 2.7 km (Fig. 20). Similarly, we define the average inflow velocity into each circle, *U*_{in}(*ρ*), as (−*dA*/*d*t)_{ob}/*πρ*. The average inflow velocity increases almost linearly with *ρ* from negative values (signifying outflow) in the inner cell to strong mean inflow winds of 16 m s^{−1} at 5 km (Fig. 21). The average wind speed $\u2061(Uin2+Vtang2)0.5$ ranges from 14 to 18 m s^{−1} as *ρ* varies from 1 to 5 km (Fig. 22). The inflow angle tan^{−1}(*V*_{tang}/*U*_{in}) decreases linearly from 120° at *ρ* = 1 km to 31° at *ρ* = 5 km (Fig. 23). Values greater than 90° signify outflow.

The corresponding circles at 2311 UTC are centered on the TVS rather than the mesocyclone and are at a moderately different height above radar level (1170 m instead of 812 m). The circulation (Fig. 18) and the average tangential velocity (Fig. 20) at *ρ* = 1 have both doubled since 2303 UTC, yet the circulation is roughly unchanged for *ρ* ≥ 3 km. It seems that angular momentum is being advected inward toward the rotation axis and vertically as well as in the Davies-Jones (2008) numerical model. The areal contraction rate (Fig. 19) and the average inflow velocity (Fig. 21) have increased by over 4 times at *ρ* = 3 km and by 50% at *ρ* = 5 km. The slight divergent flow at small *ρ* may be associated with the vortex core spreading with height or with centrifuging of debris. The large average inflow velocities of 24 m s^{−1} at *ρ* = 4.5–5 km (Fig. 21) indicate a surge of high momentum air (Bluestein et al. 2015) toward the TVS. The average windspeeds (Fig. 22) have increased at *ρ* = 1–2 km owing to stronger tangential winds and at *ρ* = 4.5–5 km owing to greater radial inflow. The inflow angle (Fig. 23) still decreases linearly with *ρ* but is less at all *ρ* than previously. High values and large temporal increases in Doppler circulation (Fig 18) and areal contraction rate (Fig. 19) clearly reveal the twin threats of intensifying tornadoes and strong larger-scale winds.

## 7. Summary

We have derived formulas that can be used for calculating the Doppler circulation around and the contraction rate of an area within a surface of constant launch angle from the field of mean Doppler velocities. Standardized measures of vortex circulation and areal contraction rates are found by bilinearly interpolating data to chosen points on within-surface circles of fixed radii concentric with circulation centers. Since the Doppler component of target velocity is the only one observed, we can compute these quantities only partially. For axisymmetric flows at ranges much greater than the circle diameter the observed circulation around and contraction rate of a horizontal circle concentric with the flow are one-half the actual values. Thus, doubling the observed values obtained from formulas herein provides useful estimates of these quantities.

We computed the Doppler circulation and mean convergence at a height of 3.5 km and radial distances of 1–3 km from the 24 May 1973 Union City TVS and found them to be practically constant. After doubling to compensate for the unobserved wind component being set to zero, the circulation value at 3.5 km AGL (~1.1 × 10^{5} m^{2} s^{−1}) agreed closely with the photogrammetrically measured value near the ground. The results indicate that the mature tornado was embedded in a region of nearly uniform strong convergence (~5.5 × 10^{−3} s^{−1} after doubling) that was about 6 km in diameter without a mesocyclone present at the time. We reproduced these results using a simple model of a convergent vortex as input to a Doppler radar emulator. We then moved the model vortex to show that for a WSR-88D with superresolution the circulation around a circle of radius 2.5 (1.5) km declines by less than 20% for ranges up to 225 (135) km.

We also analyzed data collected by the KOUN WSR-88D on the 31 May 2013 El Reno storm. At the time of tornado formation, the analysis reveals a two-celled mesocyclone with strong inflow at 5 km from the circulation axis. Eight minutes later the situation had worsened considerably. Large intensifications in circulation near the axis and in the influx of momentum at the 5 km radius signified greater probability of violent tornadoes surrounded by strong mesocyclonic and storm-scale winds.

Work still needs to be done to statistically test the usefulness of circulation and areal expansion rate in tornado warnings. We plan on making the observed circulations more accurate by using Xu et al.’s (2017) method to better position the vortex axis within an ASI. We will also evaluate the effect of noise on the simulated circulation measurements.

## Acknowledgments

Drs. Qin Xu and Jeffrey Snyder of the National Severe Storms Laboratory (NSSL) and the three formal reviewers provided valuable comments. Thanks also to Mr. Charles Kuster of NSSL for assisting with the dealiasing of the 31 May 2013 KOUN WSR-88D data.

### APPENDIX

#### Explanations for Doubling Observed Circulation and for Quadrupole Patterns

Here we show why the actual circulation and the areal contraction rate are about twice the Doppler observed values and we also reveal the origins of the spurious quadrupole patterns in the mean Doppler velocity fields of axisymmetric flows. For these purposes horizontal plane geometry with Cartesian coordinates (*x*, *y*) and orthonormal basis vectors **i** and **j** suffices. Assume axisymmetric flow around the origin and a Doppler radar with perfect resolution at *x* = −∞, *y* = 0 so that we do not have to calculate mean Doppler velocities and to account for divergence of the radials. Define polar coordinates (*σ*, *θ*) where *x* = *σ* cos*θ* and *y* = *σ* sin*θ* and conversely

The velocity vector is

where

are the orthonormal basis vectors of the polar coordinates. Here *U* and *V* are the axisymmetric flow’s radial and tangential velocities, *u* is the Doppler velocity and *υ* is the unobserved velocity component where

The Doppler circulation is therefore

which is one-half the actual circulation around *C*. With the same assumptions, the proof that the Doppler areal contraction rate is one-half the actual areal contraction rate is similar.

By differentiating (A4) using the chain rule, differentiating (A1) to obtain ∂*σ*/∂*x*, ∂*σ*/∂*x*, ∂*θ*/∂*x*, and ∂*θ*/∂*y*, and using trigonometric identities, we obtain for axisymmetric flow

where here ∂*u*/∂*x* is the Doppler (observed) divergence, ∂*υ*/∂*y* is the unobserved divergence, ∂*υ*/∂*x* is the unobserved vorticity, and −∂*u*/∂*y* is the Doppler vorticity. The circulation around (areal expansion rate) of a cell is the mean cell Doppler vorticity (divergence) times the cell area, which is practically constant near a flow axis far from the radar. The terms in 2*θ* in (A12) and (A9) give rise to spurious quadrupole patterns in the fields of cell circulations and cell expansion rates, respectively. Only in a flow *U* = −*cσ*, *V* = Ω*σ* in solid-body rotation with constant convergence 2*c* and vorticity 2Ω are the quadrupoles entirely absent.

## REFERENCES

*Wea. Forecasting*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*J. Atmos. Oceanic Technol.*

*Proc. Symp. on Tornadoes: Assessment of Knowledge and Implications for Man*, Lubbock, TX, Texas Tech University, 151–174

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*28th Conf. on Radar Meteorology*, Austin, TX, Amer. Meteor. Soc., 14–16

*J. Atmos. Oceanic Technol.*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*Wea. Forecasting*

*Doppler Radar and Weather Observations*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*Wea. Forecasting*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*Wea. Forecasting*

*The Mathematics of Physics and Chemistry*

*Wea. Forecasting*

*Motion and Motion Systems*. Vol. I,

*Weather Analysis and Forecasting*, McGraw-Hill, 428 pp

*Z. Angew. Math. Phys.*

*Bull. Amer. Meteor. Soc.*

*Wea. Forecasting*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*Wea. Forecasting*

*Bull. Amer. Meteor. Soc.*

*J. Atmos. Oceanic Technol.*

## Footnotes

^{a}

Emeritus.

^{1}

The areal contraction rate is the negative of the areal expansion rate.

^{2}

Except at exceedingly short range, cell areas near a vortex are locally constant.

The Tornado: Its Structure, Dynamics, Prediction, and Hazards,Geophys. Monogr., Vol. 79, Amer. Geophys. Union, 203–211