J Atmos Ocean Technol 33(12): 2615-2638

Making Aircraft Vortices Visible to Radar by Spraying Water into the Wake

1. Introduction

a. Motivation

Aircraft trailing vortices are a hazard to following airplanes during take-off and landing (e.g., Barbagallo 2014). For a review of their dynamics, see Spalart (1998). Due to their mutually induced velocity, a pair of trailing vortices generally descends from the altitude where it was generated. (For exceptions to this in strongly stratified conditions, see Spalart, 1996). Therefore, a wake encounter may occur when a following aircraft finds itself below¹ the path of the leading aircraft. The possibility for this is increased during take-off and landing, and the present work addresses the latter situation. When cleared for a visual approach to landing, the pilot of the following aircraft can visually attempt to remain above the path of the leader by, for example, flying at a higher glide slope to the same touchdown point as the leader. Even in visual approaches, however, things can and do go awry (Barbagallo 2014, §2.3). This means that efforts to find an all-weather wake sensor should be continued.

When the ceiling is less than 1000 ft and the visibility less than 3 statute miles, operations must be conducted under Instrument Flight Regulations (IFR). In this case, air traffic controllers maintain separations according to the weight categories of the leading and following aircraft; see Table 1. These separations have started to limit capacity at some airports (Crouch et al. 2001b) and we refer the reader to a report (Broderick et al. 2008) by a committee of the National Research Council entitled “Wake Turbulence: An Obstacle to Increased Air Traffic Capacity.” About two decades ago, NASA initiated the AVOSS (Aircraft Vortex Spacing System) program whose aim was to make aircraft spacing dynamic through a combination of vortex sensing and real-time flowfield simulation. The present work is motivated by the vortex detection aspect of the AVOSS program. A current effort that has similar aims as the defunct AVOSS program is WakeNet-Europe (www.wakenet.eu). Finally, we mention the development of a system to hasten the breakup of trailing vortices by exciting vortex instabilities through periodic motion of control surfaces (Crouch et al. 2001a,b). Such a system could be implemented together with one for wake detection.

TABLE 1

IFR separation standards (in nautical miles) for arrivals on the same runway (Barbagallo 2014).

	Follower
Leader	Super	Heavy	B757	Large	Small
Super	3	6	7	7	8
Heavy	3	4	5	5	6
B757	3	4	4	4	5
Large	3	3	3	3	4
Small	3	3	3	3	3

b. Previous work using radar to detect aircraft wakes

The main vortex detection technology tested by the AVOSS program was infrared lidar. A concern raised about lidar was that since water vapor strongly absorbs infrared, it would not be usable in IFR conditions. Another concern is that optical systems are more expensive and difficult to maintain than radar. These concerns motivate consideration of radar.

1) Clear air reflectivity

The first observational and theoretical efforts on radar detection of wakes were in the context of clear air. Systematic tests by Gilson (1992) showed that the wake of a C-5A aircraft could be detected by radars having 2–7 MW of peak power at frequencies from 0.162 to 5.67 GHz, with no return at 35 GHz. Shariff and Wray (2002) analyzed the reflectivity in Gilson’s test using a model of a vortex pair descending in a stratified atmosphere, carrying with it an oval of atmospheric air from the altitude at which the wake formed. This leads to a gradient in refractive index between the oval and ambient air. Another mechanism investigated by Shariff and Wray (2002), which has peak reflectivity at a low frequency of 50 MHz, is the pressure (hence density) gradient in each vortex. Both of these mechanisms have some drawbacks for practical use: (a) the radar cross-section is small (−60 to −80 dB m^{); (b) the first mechanism depends on atmospheric stratification, which has seasonal, geographic, and diurnal variations; and (c) the frequency of the second mechanism is the same as that of Stratospheric-Tropospheric radars, which have an antenna array on the ground looking upward. For the present application, technology would have to be developed for aiming the radar by use of phasing.}Li et al. (2011) extended the work of Shariff and Wray (2002) with a better calculation of the compressibility-induced variation in air density in the wake. More importantly, they showed that the atmospheric gradient of water vapor is wound up into a spiral by each vortex, which allows scattering at high frequencies.

Babaresco et al. (2008) conducted tests using an X-band (9.6 GHz) radar for aircraft taking-off and flying straight and level at 1500 m altitude. For the take-off case, the range was about 700 m looking roughly sideways using peak pulse power of only 20 W. For the aircraft at an altitude of 1500 m, the radar was looking straight up and the peak pulse power was 75 W. A time doppler plot indicated a spiral structure in each vortex. However, in staring mode, a return is obtained for only five seconds, at most. Babaresco et al. (2008) did not report how far downstream of the aircraft the wake could be detected; it appears that this distance is very short.

In conclusion, clear air reflectivity is an interesting prospect, but more careful and better documented observational campaigns combined with theoretical efforts are needed.

2) Exploiting natural precipitation

A different approach, which the author first heard about from Robert Neece (NASA Langley, personal communication) around 2000, exploits the fact that water droplets (in the form of fog or rain) are present in IFR weather conditions. They are strong radar reflectors, and, when an airplane wake sets these particles into motion, they can be separated from ambient droplets via their doppler signature. This strategy has the advantage that standard doppler weather radars could be used. Seliga and Mead (2009) demonstrated the feasibility of the approach using a W-band (94 GHz) radar in an opportunistic test with only 100 mW of peak power. An analysis of the reflectivity of this mechanism has been conducted by Liu et al. (2013), and further measurements are reported in Babaresco (2012). The main drawback of this approach is that there may be long stretches along the flight path where a sufficient number of natural droplets is unavailable.

c. Present approach

The approach proposed here is related to the use of natural precipitation for wake detection discussed in the previous paragraph, except that, to provide a persistent radar target, water spray is injected into the wake. We envision one nozzle on each side of the airplane that injects the water spray near the wing trailing edge at a specified spanwise location. One possibility for nozzle placement is the aft tip of a flap track fairing. A pump and water tank could be located nearby within the wing.

The main constraint imposed by nature is droplet retention in the wake: a nozzle produces a distribution of droplet sizes; the smallest ones evaporate away while the largest ones fall out of the wake due to gravity. However, IFR conditions are correlated with high humidity (§2l). This reduces the rate of evaporation and makes the approach feasible. In a non-IFR case we consider that has moderate humidity, only the Ka and W band radars give sufficient reflectivity in a patch above each vortex. However, it is shown that spectral processing reduces noise and allows detection even when the signal to noise ratio in individual pulse returns is < 1. Finally, there is no reason why the present approach and that of using natural precipitation could not be combined using the same radar.

a. Motivation

TABLE 1

IFR separation standards (in nautical miles) for arrivals on the same runway (Barbagallo 2014).

	Follower
Leader	Super	Heavy	B757	Large	Small
Super	3	6	7	7	8
Heavy	3	4	5	5	6
B757	3	4	4	4	5
Large	3	3	3	3	4
Small	3	3	3	3	3

b. Previous work using radar to detect aircraft wakes

1) Clear air reflectivity

In conclusion, clear air reflectivity is an interesting prospect, but more careful and better documented observational campaigns combined with theoretical efforts are needed.

2) Exploiting natural precipitation

1) Clear air reflectivity

In conclusion, clear air reflectivity is an interesting prospect, but more careful and better documented observational campaigns combined with theoretical efforts are needed.

2) Exploiting natural precipitation

c. Present approach

2. Calculation methods

a. General procedure

We use aircraft centered coordinates (x,y,z), where x is streamwise, y is spanwise, and z is vertical. Air and water are denoted by subscripts ‘a’ and ‘w’, respectively.

The calculation has three steps. The first creates a sample of droplet radii a(t) from a size distribution pertinent to aerial spray nozzles. The second step generates a spray trail on the starboard side of the aircraft wake. This involves injecting droplets from the sample into the wake, tracking the position X(t) of each droplet and its radius a(t) as it evaporates. The flow-field model consists of two counter-rotating vortices whose height decreases with downstream distance x behind the aircraft. The third and final step mirrors the starboard trail to the port side and, for a given set of radar parameters, computes the reflectivity of both trails together.

Since it is prohibitive to track all of the actual droplets, the spray trail computed in step two consists of a certain number, N_comp, of computational droplets. In the reflectivity calculation, each computational droplet is taken to represent a multiplicity, M_true, of actual droplets, which is simply the ratio of the desired injected volume to the volume injected in the computation.

The equations for droplet motion and evaporation are integrated using the routine LSODE which is described in Radhakrishnan and Hindmarsh (1993) and available from NETLIB. The routine chooses a time step based on a specified error tolerance. Since the flow-field model we have adopted is steady in a reference frame moving with the aircraft, it is not necessary to inject droplets at every time step, which would continuously increase the number of computational droplets. Instead, droplets are injected only during a certain time interval Δt₀ at the beginning of the calculation. This set of droplets is then evolved for successive Δt₀ intervals and appended to a file at the end of every Δt₀ interval. At the end of the computation, the file contains a spray trail that is 7 nm miles long, whose radar reflectivity we subsequently analyze. In actual practice, an aircraft would likely not generate a trail of several nautical miles behind it. Rather, it would release the spray for a short period at pre-selected locations during its approach.

A detailed description of each part of the procedure is described in the following subsections. The evaporation calculation is described in Appendix A2. To avoid stiffness of the system of evolution equations, droplets are removed from the calculation when their radius becomes < 20 μm; at this point their reflectivity is too small to significantly affect the received power.

b. Droplet trajectory

The position X(t) and velocity U(t) of a droplet of mass m_d evolves according to

\frac{d X}{d t} = U (t),

(1)

\frac{d U}{d t} = F_{D} / m_{d} - g_{eff} \hat{Z},

(2)

where g_eff = (1 − ρ_a/ρ_w)g is the effective gravity accounting for buoyancy. The drag force F_D is given by

F_{D} = C_{D} \frac{1}{2} ρ_{a} {| u_{rel} |}^{2} π a^{2} \frac{u_{rel}}{| u_{rel} |},

(3)

where

u_rel = u(X) − U,

(4)

is the velocity of the air flow, u(X), relative to the droplet. Evaluation of the drag coefficient, C_D, is described in Appendix A1.

c. Droplet size distribution of aerial spray nozzles

When the water jet issues from the nozzle, it will encounter a blast of free-stream air with a speed of 77.2 m s⁻^{. A comparable situation in the literature is that of a cylindrical liquid jet surrounded by an annulus of co-flowing air (}Lorenzetto and Lefebvre 1977;Varga et al. 2003) which shows that much smaller droplets are produced than for the case of still air (with the same velocity of the liquid jet). This is due to the occurrence of both the Kelvin-Helmholtz instability and Rayleigh-Taylor instabilities. The former is driven by the shear between the air and liquid flow and leads to a smaller instability wavelength. The latter arises due to the acceleration of liquid droplets by the drag force of the air stream. If the droplets are too small, they quickly evaporate. To produce larger droplets, both instabilities can be mitigated by reducing the relative velocity between the liquid jet and free-stream air. This is accomplished by increasing the driving pressure. However, if the droplets are too large, they fall out of the wake. Hence, there is an optimum droplet size.

We were fortunate that an experimental study (Fritz and Hoffmann 2015), which uses a wind-tunnel to mimic the free stream air flow, has recently been performed to characterize the droplet size distribution produced by nozzles used for aerial agricultural spraying. This study did indeed show that larger pressures produce larger droplets. Dr. B. Fritz kindly provided us with an Excel program, developed from that study, which gives parameters of the droplet size distribution for various nozzles, free-stream speeds, and driving pressures. Use of these parameters is now described.

Let p(a) be the probability density such that p(a)da is the probability that the droplet radius is in the interval [a, a + da]. The log-normal distribution

p (a) = \frac{1}{\sqrt{2 π} a σ} e^{- {ln}^{2} (a / a_{0}) / 2 σ^{2}},

(5)

with parameters a₀ and σ, is commonly used in the spray literature. Fritz’s Excel program provides information about the function Q(a), defined to be the fractional volume occupied by droplets of radius ≤ a. One can show from appropriately integrating (5) that

Q(a) = ½(1 + erfξ),

(6)

where

ξ \equiv \frac{1}{\sqrt{2} σ} [ln (a / a_{0}) - 3 σ^{2}] .

(7)

The Excel program provides a₀_.₅ and a₀_.₉ defined such that Q(a₀_.₅) = 0.5 and Q(a₀_.₉) = 0.9. Using them, (6) can be numerically inverted to yield the parameters a₀ and σ of the log-normal distribution (5). The Excel program also provides a₀_.₁, which we did not use because we wished to nail the size distribution for large droplets which contribute most to reflectivity.

d. Choice of nozzle

What is the best drop size distribution? A set of N_d droplets over which the incident beam is assumed to have uniform intensity, has a reflectivity proportional to (e.g., Doviak and Zrnić 1984, p. 58)

ζ = \sum_{i = 1}^{N_{d}} a_{i}^{6},

(8)

assuming Rayleigh scattering. Maximizing this subject to fixed volume of water and fixed N_d gives the result that all droplets must be of the same size. Given this result, maximizing ζ with respect to the number N_d subject to fixed volume gives N_d = 1, i.e., all the volume must be in one droplet. However, such a droplet would likely fall too rapidly. To minimize droplet loss by sedimentation, droplets must not have a terminal velocity larger than the vortex descent speed, W_descent = 1.75 m s⁻^{in the present case. Consulting the terminal velocity plot in}Pruppacher and Klett (1997, p. 416) we conclude that a must be ≤ 200 μm. Hence the best distribution is uniform with a drop radius a = 200 μm. Since the rate of droplet evaporation is ∝ 1/a, i.e., small droplets evaporate faster than larger ones, the above conclusion is not altered by including evaporation.

The above considerations suggest that the following “rate of ζ” could be used to initially evaluate different nozzles without having to perform an evaporation and reflectivity calculation:

{\dot{ζ}}_{nozzle}^{<} \equiv \frac{1}{Δ t} \sum_{i : a_{i} < 200 μ m}^{N_{d} (Δ t)} a_{i}^{6},

(9)

where the < superscript denotes the exclusion from the sum of droplets that fall away and N_d(Δt) denotes the number of droplets produced by the nozzle in a period Δt. This suggestion will be tested in §3h.

Table 2 lists the parameters for four candidate nozzles and operating conditions that were selected from B. Fritz’s Excel program based on having peak probability density at a large radius (which rarely exceeded 100 μm) and a high flow-rate. The flow-rate for nozzle 1 was provided by Calvin Kroes (private communication) of CP Products. For nozzles 2 and 3, flow-rates were extrapolated from the values of 5.30 gpm and 2.45 gpm, respectively, at 60 psi reported on the manufacturer’s data sheets (www.translandllc.com/wp-content/uploads/2015/08/Aerial-Flow-Chart-20152.pdf and www.cpproductsinc.com/images/stories/downloads/Misc-Tables/A1-Web%20Aerial%20Tip%20Rate%20Chart.pdf). The extrapolations assumed a square-root dependence of flow-rate on pressure (Lefebvre 1989, p. 157) expected from Bernoulli’s principle, and which the manufacturer’s data follows well. For nozzle 4, the flow-rate versus pressure provided on the manufacturer’s website has a linear rather than square-root dependence. We are grateful to Dr. Brad Fritz (USDA) for measuring the actual flow-rate for us. It turned out to be much lower than the value provided by on the website.

TABLE 2

Aerial nozzle parameters for the operating conditions specified. Notes: 1: ‘Str’ denotes a straight-stream. 2: Values obtained from B. Fritz’s Excel program. 3: See text for how flow-rates were obtained.

	Nozzle 1	Nozzle 2	Nozzle 3	Nozzle 4

Model	CP-09	CP-09	CP11TT	Davidon-Triset
Pressure (psi)	90	90	90	90
Airspeed (mph)	175	175	175	175
Deflection-plane/body angle	0°	0°	0°	0°
Fan angle	Str.¹	Str.	Str.	Str.
Orifice Code			20
Orifice diameter (in)	0.125	0.172	0.105	0.125
a_{0 5} (μm)^{Note 2}	178.5	146.5	183	239.5
a_{0 9} (μm)^{Note 2}	315	303	359	463
a₀ of log-normal	99.02	55.83	79.85	108.3
σ of log-normal	0.443	0.567	0.526	0.514
Flow-rate (gpm)^{Note 3}	3.70	6.49	3.00	3.06
U_exit (m s⁻⁾	29.5	27.3	33.9	24.4
No. of nozzles per side	1	1	1	1
Gallons per nm (two sides)	2.96	5.19	2.4	2.45

e. Droplet injection in the computation

Let the origin of coordinates be in the symmetry plane of the aircraft (which corresponds to y = 0) with same axial and vertical location as the droplet injector. Droplets are injected in a grid pattern within a square of width w_square in the yz-plane. The pattern consists of n_square × n_square droplets. The square is centered at (x,y,z) = (0, fb/2,0), where f represents the fractional spanwise distance from the aircraft center plane to the wing-tip, and was chosen to be f = 0.5. A square shape was chosen because the nozzles we have selected are not of the flat fan type. The streamwise extent of the computed mist trail was chosen to be ℓ_trail = 7 nm (168 seconds of elapsed time from injection) since we wish to detect the trail at 6 nm, the longest distance for which one would want to detect the wake of a heavy aircraft under current separation rules; see Table 1.

The square pattern is injected n_x times and the time interval between injections is Δt_inject, whose value is is chosen so that the x spacing between droplets is the same as in the cross-sectional (yz) plane. After a time period Δt₀ ≡ n_xΔt_inject, an n²_square ×n_x slab of particles has been injected, which is then advanced for successive Δt₀ periods to form the entire trail. Inertial particles with a small Stokes number tend to an attractor (Haller and Sapsis 2008;Sapsis and Haller 2010) independent of injection location, and therefore where droplets end-up should be insensitive to where they are injected. A brief check on insensitivity to initial conditions will be presented in §3b.

The initial velocity of droplets is set equal to the air velocity, which is justified as follows. From the equations of droplet motion, (2)–(4), the characteristic relaxation time for a droplet to start following a new air speed, imposed at t = 0, say, is

τ_{relax} \equiv {| \frac{1}{u_{rel}} \frac{d u_{rel}}{d t} |}_{0}^{- 1} = \frac{16}{3} \frac{ρ_{w}}{ρ_{a}} \frac{a^{2}}{ν_{a}} {(Re C_{D})}_{0}^{- 1} .

(10)

where the subscript ‘0’ means that the quantity is calculated at t = 0. From this, the characteristic relaxation distance, ℓ_relax = u_rel(0)τ_relax can be evaluated. Note that the initial air speed relative to the water jet is given by u_rel(0) = U_app −U_exit, where U_app is the approach velocity of the aircraft given in Table 3, and U_exit is the exit velocity of the water jet given in Table 2. The experiment of Fritz and Hoffmann (2015) measured the size distribution 1.8 m downstream of the nozzle for all straight-stream nozzles. This value is from a private communication from B. Fritz and represents a correction from a value of 1.5 m reported in Fritz and Hoffmann (2015). Figure 1b plots ℓ_relax as a function of drop radius for nozzle 1. Inspecting it together with the size distribution in Figure 1a, we conclude that most of the droplets are following the air stream at the measurement station of the experiment. If this had not been the case and there had been a relative velocity large enough to give Weber numbers ≳ 10, then it would have been necessary to model further droplet break-up using a secondary break-up model (e.g., Apte et al. 2003).

An external file that holds a picture, illustration, etc.
Object name is nihms864032f1.jpg

Fig. 1

(a) Probability density function (pdf) of the drop sizes produced by nozzle 1 with the conditions listed in Table 2. The result is based on parameters provided by B. Fritz’s Excel program. (b) Distance required for a droplet of a given size to begin moving with an imposed air flow.

TABLE 3

Parameters of a typical heavy aircraft.

Parameter	Value
Weight, W	500,000 lb
Wing span, b	60 m
Vortex spacing, b₀	47.9 m
Vortex circulation, Γ	526 m^s⁻¹
Approach speed, U_app	150 knots

f. Flowfield of two counter-rotating wing-tip vortices

We consider an aircraft flying straight and level at altitude z₀ = 0 in aircraft coordinates. The velocity field of the airflow in the wake is denoted by lower case u(x). This velocity field consists of the free stream, $U_{app} \hat{x}$ (where U_app is the approach speed of the aircraft), superposed with the flow induced by a pair of counter-rotating vortices with circulations ±Γ. The centerline of each vortex is at spanwise location y_vort = ±b₀. Due do their mutually induced velocity, the height z_vort(x) of the vortex pair decreases with distance x behind the wing as follows:

z_vort(x) = z₀ − W_desc^t,

(11)

where W_desc = Γ/2πb₀ is the descent speed of the vortex pair. The quantity t = x/U_app is time since the vortex at x was shed from the wing.

Each vortex induces a circumferential velocity u_θ(r) in the cross-plane (yz). For u_θ(r), a profile fit to flight data by P. Spalart (Private communication) of Boeing is used:

u_{θ} = \frac{Γ}{2 π r} {\begin{cases} 1188.59 η^{2}, & η < 0.0103; \\ {[1 + {(1.27 + 0.25 log η)}^{- 14}]}^{- 1 / 14}, & otherwise; \end{cases}

(12)

where η ≡ r/b₀.

For an elliptically loaded wing, lifting line theory (Batchelor 1967) gives the vortex spacing as

b_{0} = \frac{π}{4} b,

(13)

where b is the wingspan, and the vortex circulation as

Γ = \frac{W}{ρ_{a} U_{app} b_{0}},

(14)

where W is the aircraft weight. We use the parameters for a typical heavy aircraft given in Table 3.

g. Calculation of the received signal and power

A pulse-doppler radar transmits a train of square wave pulses that modulate a carrier wave of frequency f = 2π/ω. The duration of each pulse is τ and the pulse repetition frequency is f_PRF. After each pulse is transmitted, the transmit-receive switch is set to the receive position and the incoming signal is sampled. Each sample at time t is said to come from the range gate r = c(t − t_t)/2, where t_t is the transmit time of the pulse and the factor of two accounts for the round-trip. Throughout, we consider the case where the transmitting and receiving antenna are the same, the so-called mono-static case.

1) Resolution shell

An important concept is that of the resolution volume, ℛ_τ(t) R_τ(t), at time t, associated with a single pulse of finite duration, τ (e.g., Yuter 2003, p. 1836). It is defined as the volume from which a signal is received at a fixed time t due to scattering by the pulse. Let t_t mark the beginning of the pulse at the transmitter. A signal from a scatterer at distance r will be received in the time interval

t − t_t = 2r/c + ξτ, 0 ≤ ξ ≤ 1,

(15)

where ξ = 0 corresponds to leading edge of the pulse and ξ = 1 to its trailing edge. Solving (15) for r gives

r = c(t − t_t − ξτ)/2, 0 ≤ ξ ≤ 1.

(16)

Equation (16) defines a spherical shell (called the resolution volume) from which a signal is received at the fixed time t. The next subsection describes how the complex voltage received at a given time is evaluated by summing the complex voltages from each droplet in the resolution volume.

2) Received power and signal-to-noise ratio

The material in this subsection is adapted from the texts Doviak and Zrnić (1984) and Ishimaru (1978). Let the “voltage” of the transmitted pulse be (the real part of)

V_{t} (t) = {\begin{matrix} V_{0} e^{i ω t}, & t \leq τ; \\ 0, & othrerwise; \end{matrix}

(17)

where V₀ is a complex amplitude, τ is the pulse width, and “voltage” is defined such that the instantaneous transmitted power is P_t(t) = V_t(t)V_t^∗(t). The voltage at the input terminals of the receiver is (the real part of) the following summation over droplets:

V (t) = \sum_{m = 1}^{N_{d}} A_{m} (t) e^{i ω_{m}} (t - 2 r_{m} / c),

(18)

where A_m is a complex scattering amplitude and

ω_m ≡ ω(1 − 2u_m/c)

(19)

is the twice doppler-shifted frequency (u_m being the radial velocity of the mth droplet) and r_m is the distance to each droplet. The summation in (18) is taken over the N_d droplets in the resolution volume associated with the pulse.

The complex scattering amplitude due to each droplet is

A_{m} = {[\frac{λ^{2} ℓ_{w} ℓ_{B}}{{(4 π)}^{3}} \frac{G^{2} (θ_{m})}{r_{m}^{4}} σ_{b m}]}^{1 / 2} V_{0} exp (i ϕ_{m}) .

(20)

The first factor in (20), namely [.]^{, is copied from the square root of the radar equation (e.g.,}Doviak and Zrnić 1984, p. 34) for power, where λ is the wavelength, and σ_b_m is the back-scattering cross-section of each droplet. The function G(θ_m) is the gain function of the antenna at the droplet, which we have assumed to depend on its angle θ_m from the beam centerline. Two loss factors (< 1) have been included in (20): ℓ_w is the two-way waveguide loss and ℓ_B is the loss due to finite bandwidth of the receiver; values assumed for the present study are given in Table 4. The last factor in (20), exp(iϕ_m), accounts for the phase shift induced by back-scattering at the mth droplet. The back-scattering cross-section and phase-shift (σ_b_m and φ_m, respectively) will be discussed further in §2k.

TABLE 4

Radar parameters. EEC: Enterprise Electronics Corp; PRF: Pulse repetition frequency. Note 1: nnn is the power in kW (250, 350, 500, or 1000). Note 2: The only value provided to us was for MIRA-35. The value for the other radars was assumed to be the same. Note 3: The value for a perfectly matched filter has been assumed for all radars. Note 4: For a rainfall rate of 12.5 mm/hr (medium to heavy rain). Note 5: MIRA-35 has full sensitivity beyond a range of 360 m (Matthias Bauer-Pfundstein, Private communication).

Manufacturer	EEC	EEC	EEC	Metek	ProSensing
Series	DWSR	DWSR	DWSR
Model	8501S	nnn1C^{Note 1}	2001X	MIRA-35	W-SACR
Frequency, f (GHz)	3	5.9	9.6	35.1	93.9
Peak power, P_t (kW)	850	250–1000	200	30	1.7
Reflector diameter (m)	4.2	4.2	2.4	1.2,2.0	0.9
½-power beam width, θ_b	1.83°	0.95°	0.95°	0.52°,0.31°	0.30°
Antenna gain, G (dB)	39.5	45	45	50.4,53.5	54.5
Pulse width, τ (μs)	0.4–2	0.2–3	0.2–2	0.1,0.2,0.4	0.05–2
Range resolution, cτ/2 (m)	60–300	30–450	30–300	15,30,60	7.5–300
PRF (kHz)	0.2–2.4	0.2–2.4	0.2–2.4	2.5,5,10	≤ 20
U_max = c PRF/4 f (m s⁻⁾	5–60	2.5–31	1.6–19	5.3,11,21	≤ 16
Receiver noise figure (dB)	2.0	2.0	2.0	6.2	6.0
2-way waveguide loss^(dB)	0.8	0.8	0.8	0.8	0.8
Finite bandwidth loss^(dB)	1.8	1.8	1.8	1.8	1.8
Rain attenuation^(dB/km)	0.005	0.03	0.12	2	7
Minimum range (m)				150^{Note 5}

Radar receivers have electronics that can obtain the real and imaginary parts (denoted I and Q) of V(t) and compute the instantaneous received power

P_r(t) = I(t)⁺Q(t)⁼V(t)V⁽t).

(21)

We note in passing that I(t) and Q(t) are the components of the real part of V(t) that are in-phase and 90° out-of-phase with the transmitted carrier, respectively. Substituting (18) into (21) and splitting the sum into two parts, following Doviak and Zrnić (1984, §4,1), gives

P_{r} (t) = \sum_{m, n} A_{m} A_{n}^{*} exp [i (ω_{m} - ω_{n}) t] exp | [- 2 i (k_{m} r_{m} - k_{n} r_{n})],

(22)

= \sum_{m} A_{m} A_{m}^{*} + \sum_{m, n, m \neq n} A_{m} A_{n} exp [i (ω_{m} - ω_{n}) t] exp [- 2 i (k_{m} r_{m} - k_{n} r_{n})],

(23)

where k_m ≡ ω_m/c. Arguments for using only the first term in (23) in order to evaluate reflectivity, i.e., for summing the powers received from each droplet, are given by Rayleigh (1945, p. 37), Beckmann (1962), and Doviak and Zrnić (1984, §4,1). The important point, which was phrased eloquently by Rayleigh, is that it is not correct to say that the power in a single return from a random distribution of droplets is the sum of the powers scattered by each. Rather, the result is true only when a large ensemble of returns from a statistically stationary target are averaged. This is most easily seen when we consider the case when all the A_m are equal (to unity, say). Then, the first term T₁ in (23) is T₁ = N_d. If droplet distances are randomly distributed in the resolution shell (assumed to be several wavelengths wide), then the magnitude of the second term will be the average of the summands times the number of terms, i.e., $T_{2} \approx N_{t}^{- 1 / 2} \times N_{t} = N_{t}^{1 / 2} \approx N_{d}$ , where N_t = N_d(N_d −1) is the number of terms in the double sum. Hence both terms in (23) are of similar magnitude. To make the second term smaller than the first, an ensemble average must be taken over many pulses. If the ensemble has N_s phase-uncorrelated samples, then the second term will be $N_{S}^{1 / 2}$

smaller than the first.

In the present application, the wake descends through the beam and so the target is not stationary, strictly speaking. In §3c we will explicitly verify that, for our case, an average of over a certain number of pulses does indeed yield the first term in (23). Note that a radar set does not have direct access to the first term in (23); only we as simulators do.

To provide a measure of detectability, we will present the signal-to-noise ratio

SNR1 ≡ P_r1/P_noise,

(24)

where P_r1 is the first term in (23) and the average noise power is

P_noise = k_BT₀F_N/τ,

(25)

where k_B = 1.381 × 10⁻^{J K}⁻^{is Boltzmann’s constant,}T₀ = 290 K is a reference temperature set by convention, and F_N is the overall noise figure of the chain of components in the receiving cascade. Values for F_N and τ are listed in Table 4 for each radar set.

Finally, since each computational droplet represents a multiplicity M_actual of actual droplets, we have

\sum_{m = 1}^{N_{d}} \to M_{actual} \sum_{m = 1}^{N_{comp}},

(26)

where N_comp is the number of computational droplets in the resolution volume. Note that all M_actual copies of each computational droplet are assumed to be at the same location and therefore their scattered voltages at the receiver add constructively. This assumption does not bias SNR1 since its calculation involves summing individual scattered powers anyway. We claim that the statistics of individual pulse returns are also not affected by this assumption; this will be verified (§3d) in a computational test where the number of computational droplets is increased.

h. Calculation of the doppler spectrum

Radars calculate a doppler spectrum for a given spatial observation location by performing a Fast Fourier Transform (FFT) of complex voltage returns (at the same range gate) from a sequence of pulses separated by Δt_pulses = 1/PRF, where PRF is the pulse repetition frequency. We shall do the same for the simulations. From a series of pulse returns, V_n, n = 0,1,…N_FFT − 1, the normalized transform

\hat{V} (k) \equiv \frac{1}{N_{FFT}} \sum_{n = 0}^{N_{FFT} - 1} V_{n} e^{- i 2 π k n / N_{FFT}}, k = 0, \dots N_{FFT} - 1,

(27)

and then the power spectrum $S (k) \equiv \hat{V} (k) \hat{V} * (k)$ is computed. Note that the frequency index k corresponds to an actual frequency

k_{actual} = {\begin{cases} k, & k \leq N_{FFT} / 2; \\ k - N_{FFT}, & N_{FFT} / 2 < k \leq N_{FFT} - 1. \end{cases}

(28)

The frequency k_actual is in units of (the period of the sequence)⁻^{= (}N_FFTΔt_pulses)⁻^{, so in units of s}⁻

f(k_actual) = PRFk_actual/N_FFT.

(29)

Equating this to −2(u_doppler/c) f gives the doppler velocity associated with each k_actual. Finally, we state that we use the Hamming window (e.g., Harris 1978).

i. Antenna gain function

We assume a Gaussian beam with transmitted power flux (power per unit area) vector

S_{t} = A exp (- θ^{2} / θ_{0}^{2}) \hat{r},

(30)

where A is a coefficient, θ is the angle from the beam centerline, and

θ_{0} = {(2 \sqrt{In 2})}^{- 1} θ_{b}

(31)

in terms of the half-power full-width, θ_b. The total power crossing a sphere of radius r is

P_{t} = 4 π r^{2} \int_{0}^{π} S_{t} \cdot \hat{r} sin θ, d θ = 2 π r^{2} θ_{0}^{2} A,

(32)

for a narrow beam. Using the definition of the gain function, G(θ), we obtain

G (θ) \equiv \frac{4 π r^{2} S_{t}}{P_{t}} = \frac{2}{θ_{0}^{2}} exp (- θ^{2} / θ_{0}^{2})

(33)

j. Radars included in the study

Table 4 lists parameters of currently operational doppler weather/cloud radars considered in the present study. The S, C, and X-band radars chosen are the DWSR series manufactured by EEC (Enterprise Electronics Corporation, Enterprise, Alabama). ARC (Advanced Radar Corporation, Boulder, Co.) makes quite similar C and X-band radars, while Baron Services (Huntsville, Al.) makes similar S, C, and X-band radars. The power and beamwidth values of the C-Band TDWR (Terminal Doppler Weather Radar) deployed at many US airports is subsumed by the range of values provided by the EEC C-band radar, and is therefore not included here.

A number of descriptions of Ka-band (35 GHz) cloud radars have appeared in the literature (Hamazu et al. 2003; Görsdorf et al. 2015). In the present work, we use parameters of the MIRA-35 radar manufactured by Metek (Elmshorn, Germany) which is described in Görsdorf et al. (2015). This choice was motivated by its relatively high power (30 kW). Other Ka-band weather radars, operational at the time of writing are: (i) Scanning 2 kW radars operated by U.S. Department of Energy’s Atmospheric Radiation Measurement (ARM) Climate Research Facility (Widener et al. 2012). (ii) The Copernicus 1 kW radar at Chilbolton Observatory (UK). (iii) An airborne 25 kW multi-frequency (X, Ka, and W-band) radar developed by Prosensing that is being used by NASA’s Langley Research Center for research into the detection and avoidance of super-cooled water droplets.

The science and technology of W-band (94 GHz) radars for cloud and precipitation research is reviewed in Kollias et al. (2007). For the present work we chose the W-SACR radar, which has been developed by the U.S. Department of Energy’s Atmospheric Radiation Measurement (ARM) program (Widener et al. 2012;Kollias et al. 2014).

With increasing frequency, f, reflectivity increases as f^{in Rayleigh’s formula (ignoring Mie-scattering corrections). Furthermore, the size of the antenna required to obtain the same beamwidth is reduced. The main drawback of high frequency is increased attenuation due to precipitation between the radar and target. For example, the last entry in}Table 4 gives the attenuation rate in medium rain at W-band as 7 dB/km. A compensating factor is that when there is precipitation, the ambient humidity is also very high and so there is minimal evaporation, and, if natural precipitation is present in the wake, it will also contribute to reflectivity.

k. Mie cross-section and phase-shift

Since we have can rather large droplets in the present application and frequencies up to 94 GHz, the back-scattering cross-section and phase-shift are obtained using Mie’s formula instead of Rayleigh’s approximation. We used subroutine BHMIE, available from Prof. B.T. Draine’s website at Princeton University, and checked the results using subroutine MIEV0 developed by Dr. W.J. Wiscombe (NASA Goddard).

Some understanding of notation is required to properly use these routines. Let the incident field be of unit magnitude and polarized in the 2-direction (defined to be perpendicular to the plane containing the incident and observer directions). For a spherical target, the scattered field in the far-field is also polarized in the 2-direction and is given by

E_{s2} (r, θ) = \frac{e^{i k r}}{k r} f_{22} (ϑ),

(34)

where f₂₂ is complex, k = 2π/λ, ϑ is the angle of the observer relative to the direction of propagation of the incident wave, and r is distance from the center of the sphere. The backscattering cross-section and phase shift are obtained as

σ_{b} \equiv 4 π r^{2} \frac{{| E_{s} (π) |}^{2}}{{| E_{i} |}^{2}} = \frac{4 π}{k^{2}} {| f_{22} (π) |}^{2},

(35)

ϕ = arg(f₂₂(π)).

(36)

At the start of the reflectivity calculation at a given frequency, we tabulate the ratio σ_b/σ_b_,_Rayleigh and the difference φ −φ_Rayleigh as a function of droplet radius a. Rayleigh’s formulas are (Ishimaru 1978, p. 19)

σ_b, Rayleigh = 4|K_ε|⁽ka)⁽πa^),

(37)

ϕ_Rayleigh = arg(K_ε),

(38)

where K_ε (a complex number) is given by

K_{ε} = \frac{ε - 1}{ε + 2} .

(39)

The quantity ε(f, T) is the complex dielectric constant of water; our convention of e^{^ωt} for the time dependence requires the imaginary part of ε be positive for an absorbing material. It is a function of frequency and temperature and was evaluated using the single Debye model of Liebe et al. (1991) as implemented in subroutines available from Prof. Chris O’Dell’s website at Colorado State University. Since the droplet temperature is almost the same for all drops (the spread was 4 C at most), the dielectric constant ε is evaluated at the average temperature of all the drops in the trail.

Figure 2 displays the Mie back-scattering cross-section σ_b (normalized by the Rayleigh value) and phase-shift φ as a function of droplet radius at the five frequencies considered in this work. At the largest radii in this work, a ≈ 600 μm, the error in using Rayleigh’s cross-section is about 20% at 35.1 GHz.

An external file that holds a picture, illustration, etc.
Object name is nihms864032f2.jpg

Fig. 2

Back-scattering cross-section σ_b (normalized by the Rayleigh value) and phase-shift ϕ versus droplet radius a at the frequencies considered in this work.

l. Choice of ambient temperature and humidity

Droplet evaporation calculations require specification of the ambient temperature and humidity. For guidance on appropriate choices, METARs (Meteorological Aerodrome Reports) during 2000–2014 were downloaded from

http://mesonet.agron.iastate.edu/request/download.phtml

and processed for the five busiest airports in the U.S. Figures 3a and b show the monthly-averaged temperature and relative humidity (RH), respectively, when IFR conditions prevailed. Figure 3c shows the percentage of reports that fall into the IFR category. One sees that the average RH is always above 90%. To further synthesize this data, yearly averages were taken (Table 5). Among the five airports, LAX has the highest rate of evaporation in IFR conditions on average since it is the warmest and driest on average. Our choice is the IFR average for LAX, namely, T = 15.2 C and RH = 92.7%. Looking at the monthly data for the other four airports, this appears to be a reasonable choice for them also: it is an approximate lower bound for their monthly RH and their temperature is higher only during the summer months when IFR reports are low.

An external file that holds a picture, illustration, etc.
Object name is nihms864032f3.jpg

Fig. 3

Temperature and humidity (averaged by month) when IFR conditions prevail at the five busiest U.S. airports

TABLE 5

Yearly-averaged temperature and humidity when IFR conditions prevail at the five busiest US airports.

Airport	T (C)	RH	% IFR Reports
ATL	12.5	95.5%	5.3%
LAX	15.2	92.7%	3.5%
DFW	9.7	94.5%	1.4%
ORD	5.3	94.2%	3.7%
JFK	12.1	93.9%	4.8%

Since it is expensive for flight tests to wait for IFR conditions to occur, and it is desirable to have a wake sensor that can work in a wider variety of atmosphere conditions, we will also consider a case of lower humidity and higher temperature, namely, RH = 60% and T = 20 C.

m. Radar placement with respect to the wake

Here, we choose the radar location (x_rad,y_rad,z_rad) in aircraft-centered coordinates. Based on current wake separations (Table 1), it should not be necessary to examine a wake more than 6 nm behind the aircraft. We therefore chose to present results for the reflectivity at x = 6 nm, the worst case for droplet loss by evaporation and sedimentation. The radar is also placed at x_rad = 6 nm so it can view the x = 6 nm wake cross-section at normal incidence. Next, we assume that the aircraft is at the touchdown point. At 6 nm from the touchdown point, the altitude of an aircraft flying a 3° glide slope is H = 582 m. Therefore, the vertical coordinate of the radar is z_rad = −582 m. For purposes of this study, we assume an aircraft flying straight and level at this altitude. Initial flight tests would also presumably have the aircraft fly straight and level. In this case, since the vortex pair descends at a speed of W_desc = 1.75 m s⁻^{, its axis makes a downward angle of 1.7}° relative to the horizontal. For an aircraft on a 3° glideslope, the vortex axis would therefore be 1.3° upward from the wing. This difference in the angle of the wake axis is expected to have a very small effect on reflectivity.

To place the radar laterally with respect to the wake, we imagine several parallel approaches that are monitored by the same radar. The largest separation between parallel runways is about 5000 ft (Doyle and McGee 1998). At 6 nm from touchdown, the lateral width of the ILS (Instrument Landing System) approach is 3182 ft for a standard 5° splay, and we imagine an aircraft that has strayed to the outer edge of this zone. If the radar is placed in the middle of the two furthest runways we obtain a lateral distance of 0.67 nm. In the presence of a crosswind, we imagine that the wake would be monitored for as long as it remained between the outer edges of the left and right ILS zones. In conclusion, we select (x_rad,y_rad,z_rad) = (6 nm,−0.67 nm,−582 m) relative to the aircraft.

The elevation angle of the radar beam from this location varies between 10.6° and 17.1° as the scanned range of z on the wake center plane varies between z = −350 and −200 m (see Figure 4). Since the beamwidth of the radar likely to be used is ≤ 1°, ground reflection will be small. To significantly reduce ground and structure clutter, the radar can be placed directly under the flight path. This would require a separate radar for each parallel runway. Another issue is loss of radar sensitivity at smaller ranges; for MIRA-35 this happens for r < 360 m (Matthias Bauer-Pfundstein, private communication). However, this loss is probably offset by the increase in power from the r^{factor in}Equation (20).

An external file that holds a picture, illustration, etc.
Object name is nihms864032f4.jpg

Fig. 4

Simulated SNR1 for five radars in a range-elevation scan of the x = 6 nm cross-section behind the aircraft. Panel (a) shows droplets colored by radius in μm. IFR ambient conditions: RH = 92.7%, T = 15.2 C. Nozzle 1.

a. General procedure

We use aircraft centered coordinates (x,y,z), where x is streamwise, y is spanwise, and z is vertical. Air and water are denoted by subscripts ‘a’ and ‘w’, respectively.

b. Droplet trajectory

The position X(t) and velocity U(t) of a droplet of mass m_d evolves according to

\frac{d X}{d t} = U (t),

(1)

\frac{d U}{d t} = F_{D} / m_{d} - g_{eff} \hat{Z},

(2)

where g_eff = (1 − ρ_a/ρ_w)g is the effective gravity accounting for buoyancy. The drag force F_D is given by

F_{D} = C_{D} \frac{1}{2} ρ_{a} {| u_{rel} |}^{2} π a^{2} \frac{u_{rel}}{| u_{rel} |},

(3)

where

u_rel = u(X) − U,

(4)

is the velocity of the air flow, u(X), relative to the droplet. Evaluation of the drag coefficient, C_D, is described in Appendix A1.

c. Droplet size distribution of aerial spray nozzles

Let p(a) be the probability density such that p(a)da is the probability that the droplet radius is in the interval [a, a + da]. The log-normal distribution

p (a) = \frac{1}{\sqrt{2 π} a σ} e^{- {ln}^{2} (a / a_{0}) / 2 σ^{2}},

(5)

Q(a) = ½(1 + erfξ),

(6)

where

ξ \equiv \frac{1}{\sqrt{2} σ} [ln (a / a_{0}) - 3 σ^{2}] .

(7)

d. Choice of nozzle

ζ = \sum_{i = 1}^{N_{d}} a_{i}^{6},

(8)

The above considerations suggest that the following “rate of ζ” could be used to initially evaluate different nozzles without having to perform an evaporation and reflectivity calculation:

{\dot{ζ}}_{nozzle}^{<} \equiv \frac{1}{Δ t} \sum_{i : a_{i} < 200 μ m}^{N_{d} (Δ t)} a_{i}^{6},

(9)

TABLE 2

	Nozzle 1	Nozzle 2	Nozzle 3	Nozzle 4

Model	CP-09	CP-09	CP11TT	Davidon-Triset
Pressure (psi)	90	90	90	90
Airspeed (mph)	175	175	175	175
Deflection-plane/body angle	0°	0°	0°	0°
Fan angle	Str.¹	Str.	Str.	Str.
Orifice Code			20
Orifice diameter (in)	0.125	0.172	0.105	0.125
a_{0 5} (μm)^{Note 2}	178.5	146.5	183	239.5
a_{0 9} (μm)^{Note 2}	315	303	359	463
a₀ of log-normal	99.02	55.83	79.85	108.3
σ of log-normal	0.443	0.567	0.526	0.514
Flow-rate (gpm)^{Note 3}	3.70	6.49	3.00	3.06
U_exit (m s⁻⁾	29.5	27.3	33.9	24.4
No. of nozzles per side	1	1	1	1
Gallons per nm (two sides)	2.96	5.19	2.4	2.45

e. Droplet injection in the computation

τ_{relax} \equiv {| \frac{1}{u_{rel}} \frac{d u_{rel}}{d t} |}_{0}^{- 1} = \frac{16}{3} \frac{ρ_{w}}{ρ_{a}} \frac{a^{2}}{ν_{a}} {(Re C_{D})}_{0}^{- 1} .

(10)

Fig. 1

TABLE 3

Parameters of a typical heavy aircraft.

Parameter	Value
Weight, W	500,000 lb
Wing span, b	60 m
Vortex spacing, b₀	47.9 m
Vortex circulation, Γ	526 m^s⁻¹
Approach speed, U_app	150 knots

f. Flowfield of two counter-rotating wing-tip vortices

z_vort(x) = z₀ − W_desc^t,

(11)

where W_desc = Γ/2πb₀ is the descent speed of the vortex pair. The quantity t = x/U_app is time since the vortex at x was shed from the wing.

Each vortex induces a circumferential velocity u_θ(r) in the cross-plane (yz). For u_θ(r), a profile fit to flight data by P. Spalart (Private communication) of Boeing is used:

u_{θ} = \frac{Γ}{2 π r} {\begin{cases} 1188.59 η^{2}, & η < 0.0103; \\ {[1 + {(1.27 + 0.25 log η)}^{- 14}]}^{- 1 / 14}, & otherwise; \end{cases}

(12)

where η ≡ r/b₀.

For an elliptically loaded wing, lifting line theory (Batchelor 1967) gives the vortex spacing as

b_{0} = \frac{π}{4} b,

(13)

where b is the wingspan, and the vortex circulation as

Γ = \frac{W}{ρ_{a} U_{app} b_{0}},

(14)

where W is the aircraft weight. We use the parameters for a typical heavy aircraft given in Table 3.

g. Calculation of the received signal and power

1) Resolution shell

t − t_t = 2r/c + ξτ, 0 ≤ ξ ≤ 1,

(15)

where ξ = 0 corresponds to leading edge of the pulse and ξ = 1 to its trailing edge. Solving (15) for r gives

r = c(t − t_t − ξτ)/2, 0 ≤ ξ ≤ 1.

(16)

2) Received power and signal-to-noise ratio

The material in this subsection is adapted from the texts Doviak and Zrnić (1984) and Ishimaru (1978). Let the “voltage” of the transmitted pulse be (the real part of)

V_{t} (t) = {\begin{matrix} V_{0} e^{i ω t}, & t \leq τ; \\ 0, & othrerwise; \end{matrix}

(17)

V (t) = \sum_{m = 1}^{N_{d}} A_{m} (t) e^{i ω_{m}} (t - 2 r_{m} / c),

(18)

where A_m is a complex scattering amplitude and

ω_m ≡ ω(1 − 2u_m/c)

(19)

The complex scattering amplitude due to each droplet is

A_{m} = {[\frac{λ^{2} ℓ_{w} ℓ_{B}}{{(4 π)}^{3}} \frac{G^{2} (θ_{m})}{r_{m}^{4}} σ_{b m}]}^{1 / 2} V_{0} exp (i ϕ_{m}) .

(20)

TABLE 4

Manufacturer	EEC	EEC	EEC	Metek	ProSensing
Series	DWSR	DWSR	DWSR
Model	8501S	nnn1C^{Note 1}	2001X	MIRA-35	W-SACR
Frequency, f (GHz)	3	5.9	9.6	35.1	93.9
Peak power, P_t (kW)	850	250–1000	200	30	1.7
Reflector diameter (m)	4.2	4.2	2.4	1.2,2.0	0.9
½-power beam width, θ_b	1.83°	0.95°	0.95°	0.52°,0.31°	0.30°
Antenna gain, G (dB)	39.5	45	45	50.4,53.5	54.5
Pulse width, τ (μs)	0.4–2	0.2–3	0.2–2	0.1,0.2,0.4	0.05–2
Range resolution, cτ/2 (m)	60–300	30–450	30–300	15,30,60	7.5–300
PRF (kHz)	0.2–2.4	0.2–2.4	0.2–2.4	2.5,5,10	≤ 20
U_max = c PRF/4 f (m s⁻⁾	5–60	2.5–31	1.6–19	5.3,11,21	≤ 16
Receiver noise figure (dB)	2.0	2.0	2.0	6.2	6.0
2-way waveguide loss^(dB)	0.8	0.8	0.8	0.8	0.8
Finite bandwidth loss^(dB)	1.8	1.8	1.8	1.8	1.8
Rain attenuation^(dB/km)	0.005	0.03	0.12	2	7
Minimum range (m)				150^{Note 5}

Radar receivers have electronics that can obtain the real and imaginary parts (denoted I and Q) of V(t) and compute the instantaneous received power

P_r(t) = I(t)⁺Q(t)⁼V(t)V⁽t).

(21)

P_{r} (t) = \sum_{m, n} A_{m} A_{n}^{*} exp [i (ω_{m} - ω_{n}) t] exp | [- 2 i (k_{m} r_{m} - k_{n} r_{n})],

(22)

= \sum_{m} A_{m} A_{m}^{*} + \sum_{m, n, m \neq n} A_{m} A_{n} exp [i (ω_{m} - ω_{n}) t] exp [- 2 i (k_{m} r_{m} - k_{n} r_{n})],

(23)

smaller than the first.

To provide a measure of detectability, we will present the signal-to-noise ratio

SNR1 ≡ P_r1/P_noise,

(24)

where P_r1 is the first term in (23) and the average noise power is

P_noise = k_BT₀F_N/τ,

(25)

Finally, since each computational droplet represents a multiplicity M_actual of actual droplets, we have

\sum_{m = 1}^{N_{d}} \to M_{actual} \sum_{m = 1}^{N_{comp}},

(26)

1) Resolution shell

t − t_t = 2r/c + ξτ, 0 ≤ ξ ≤ 1,

(15)

where ξ = 0 corresponds to leading edge of the pulse and ξ = 1 to its trailing edge. Solving (15) for r gives

r = c(t − t_t − ξτ)/2, 0 ≤ ξ ≤ 1.

(16)

2) Received power and signal-to-noise ratio

The material in this subsection is adapted from the texts Doviak and Zrnić (1984) and Ishimaru (1978). Let the “voltage” of the transmitted pulse be (the real part of)

V_{t} (t) = {\begin{matrix} V_{0} e^{i ω t}, & t \leq τ; \\ 0, & othrerwise; \end{matrix}

(17)

V (t) = \sum_{m = 1}^{N_{d}} A_{m} (t) e^{i ω_{m}} (t - 2 r_{m} / c),

(18)

where A_m is a complex scattering amplitude and

ω_m ≡ ω(1 − 2u_m/c)

(19)

The complex scattering amplitude due to each droplet is

A_{m} = {[\frac{λ^{2} ℓ_{w} ℓ_{B}}{{(4 π)}^{3}} \frac{G^{2} (θ_{m})}{r_{m}^{4}} σ_{b m}]}^{1 / 2} V_{0} exp (i ϕ_{m}) .

(20)

TABLE 4

Manufacturer	EEC	EEC	EEC	Metek	ProSensing
Series	DWSR	DWSR	DWSR
Model	8501S	nnn1C^{Note 1}	2001X	MIRA-35	W-SACR
Frequency, f (GHz)	3	5.9	9.6	35.1	93.9
Peak power, P_t (kW)	850	250–1000	200	30	1.7
Reflector diameter (m)	4.2	4.2	2.4	1.2,2.0	0.9
½-power beam width, θ_b	1.83°	0.95°	0.95°	0.52°,0.31°	0.30°
Antenna gain, G (dB)	39.5	45	45	50.4,53.5	54.5
Pulse width, τ (μs)	0.4–2	0.2–3	0.2–2	0.1,0.2,0.4	0.05–2
Range resolution, cτ/2 (m)	60–300	30–450	30–300	15,30,60	7.5–300
PRF (kHz)	0.2–2.4	0.2–2.4	0.2–2.4	2.5,5,10	≤ 20
U_max = c PRF/4 f (m s⁻⁾	5–60	2.5–31	1.6–19	5.3,11,21	≤ 16
Receiver noise figure (dB)	2.0	2.0	2.0	6.2	6.0
2-way waveguide loss^(dB)	0.8	0.8	0.8	0.8	0.8
Finite bandwidth loss^(dB)	1.8	1.8	1.8	1.8	1.8
Rain attenuation^(dB/km)	0.005	0.03	0.12	2	7
Minimum range (m)				150^{Note 5}

Radar receivers have electronics that can obtain the real and imaginary parts (denoted I and Q) of V(t) and compute the instantaneous received power

P_r(t) = I(t)⁺Q(t)⁼V(t)V⁽t).

(21)

P_{r} (t) = \sum_{m, n} A_{m} A_{n}^{*} exp [i (ω_{m} - ω_{n}) t] exp | [- 2 i (k_{m} r_{m} - k_{n} r_{n})],

(22)

= \sum_{m} A_{m} A_{m}^{*} + \sum_{m, n, m \neq n} A_{m} A_{n} exp [i (ω_{m} - ω_{n}) t] exp [- 2 i (k_{m} r_{m} - k_{n} r_{n})],

(23)

smaller than the first.

To provide a measure of detectability, we will present the signal-to-noise ratio

SNR1 ≡ P_r1/P_noise,

(24)

where P_r1 is the first term in (23) and the average noise power is

P_noise = k_BT₀F_N/τ,

(25)

Finally, since each computational droplet represents a multiplicity M_actual of actual droplets, we have

\sum_{m = 1}^{N_{d}} \to M_{actual} \sum_{m = 1}^{N_{comp}},

(26)

h. Calculation of the doppler spectrum

\hat{V} (k) \equiv \frac{1}{N_{FFT}} \sum_{n = 0}^{N_{FFT} - 1} V_{n} e^{- i 2 π k n / N_{FFT}}, k = 0, \dots N_{FFT} - 1,

(27)

and then the power spectrum $S (k) \equiv \hat{V} (k) \hat{V} * (k)$ is computed. Note that the frequency index k corresponds to an actual frequency

k_{actual} = {\begin{cases} k, & k \leq N_{FFT} / 2; \\ k - N_{FFT}, & N_{FFT} / 2 < k \leq N_{FFT} - 1. \end{cases}

(28)

The frequency k_actual is in units of (the period of the sequence)⁻^{= (}N_FFTΔt_pulses)⁻^{, so in units of s}⁻

f(k_actual) = PRFk_actual/N_FFT.

(29)

Equating this to −2(u_doppler/c) f gives the doppler velocity associated with each k_actual. Finally, we state that we use the Hamming window (e.g., Harris 1978).

i. Antenna gain function

We assume a Gaussian beam with transmitted power flux (power per unit area) vector

S_{t} = A exp (- θ^{2} / θ_{0}^{2}) \hat{r},

(30)

where A is a coefficient, θ is the angle from the beam centerline, and

θ_{0} = {(2 \sqrt{In 2})}^{- 1} θ_{b}

(31)

in terms of the half-power full-width, θ_b. The total power crossing a sphere of radius r is

P_{t} = 4 π r^{2} \int_{0}^{π} S_{t} \cdot \hat{r} sin θ, d θ = 2 π r^{2} θ_{0}^{2} A,

(32)

for a narrow beam. Using the definition of the gain function, G(θ), we obtain

G (θ) \equiv \frac{4 π r^{2} S_{t}}{P_{t}} = \frac{2}{θ_{0}^{2}} exp (- θ^{2} / θ_{0}^{2})

(33)

j. Radars included in the study

k. Mie cross-section and phase-shift

E_{s2} (r, θ) = \frac{e^{i k r}}{k r} f_{22} (ϑ),

(34)

σ_{b} \equiv 4 π r^{2} \frac{{| E_{s} (π) |}^{2}}{{| E_{i} |}^{2}} = \frac{4 π}{k^{2}} {| f_{22} (π) |}^{2},

(35)

ϕ = arg(f₂₂(π)).

(36)

σ_b, Rayleigh = 4|K_ε|⁽ka)⁽πa^),

(37)

ϕ_Rayleigh = arg(K_ε),

(38)

where K_ε (a complex number) is given by

K_{ε} = \frac{ε - 1}{ε + 2} .

(39)

Fig. 2

Back-scattering cross-section σ_b (normalized by the Rayleigh value) and phase-shift ϕ versus droplet radius a at the frequencies considered in this work.

l. Choice of ambient temperature and humidity

http://mesonet.agron.iastate.edu/request/download.phtml

Fig. 3

Temperature and humidity (averaged by month) when IFR conditions prevail at the five busiest U.S. airports

TABLE 5

Yearly-averaged temperature and humidity when IFR conditions prevail at the five busiest US airports.

Airport	T (C)	RH	% IFR Reports
ATL	12.5	95.5%	5.3%
LAX	15.2	92.7%	3.5%
DFW	9.7	94.5%	1.4%
ORD	5.3	94.2%	3.7%
JFK	12.1	93.9%	4.8%

m. Radar placement with respect to the wake

Fig. 4

3. Results

a. Signal-to-noise ratio for the IFR case

We begin by considering IFR ambient conditions (RH = 92.7%,T = 15.2 C) chosen as described in §2l. Nozzle 1 from Table 2 is used and parameters for the injected square of droplets, described in §2e, are n_square = 15, n_x = 120, and w_square = 1 m. Figure 4 shows simulated values of SNR1, calculated using (24), for the five radars listed is Table 4. Each SNR1 plot is an instantaneous range-elevation scan of the x = 6 nm cross-section of the wake and each location on the plot corresponds to the mid-radius of a resolution shell along the beam centerline. The pulse width is chosen to be τ = 0.2 μs for all the radars except for DWSR-8501S, in which case the lowest available τ of 0.4 μs is used. Droplets in a 30 m thick axial slab centered at x = 6 nm are shown in panel (a). Due to centrifugation, larger droplets lie at greater distances from the vortex center, which is devoid of droplets. The very large particles sediment due to gravity after being centrifuged. Except for the S-band radar, all radars give SNR1 > 10 dB at most points surrounding the vortices; the reflectivity is higher for the higher frequency radars. The W-SACR radar gives the highest reflectivity despite having the smallest pulse power.

For all radars, there is a drop in reflectivity near the 2 o’clock and 4 o’clock positions for the left vortex (8 and 10 o’clock positions for the right vortex). This manifests as a crescent-wrench shaped reflectivity pattern that is most prominent for the DWSR 2001X radar. How this feature is related to the droplet configuration, which in turn is related to the vortex flowfield, remains to be elucidated.

Figure 5 shows that with its lowest available pulse width of 0.05 μs, the W-SACR radar is able to resolve some of the spiral structure of the droplet pattern at the expense of some loss in SNR1.

An external file that holds a picture, illustration, etc.
Object name is nihms864032f5.jpg

Fig. 5

Simulated SNR1 for W-SACR with a pulse width of 0.05 μs.

b. Insensitivity to initial condition

To test sensitivity to initial conditions, instead of injecting droplets in a regular grid pattern on each square, droplets were randomly placed in the squares. Figures 6a and b show that both the droplet configuration at x = 6 nm as well SNR1 for DWSR-2001X are changed very little; compare with Figures 4a and d. We expect this to be true for all the radars as well. In another test (Figure 6c and d), the width of the square, w_square, was reduced from 1 m to 50 cm, keeping the number of droplets fixed. This increases the initial number density in the cross-plane (yz) and, to keep the spacing the same in the streamwise (x) direction, the injection interval Δt_inject was also halved. One might think that this would increase the number density downstream. However, the flow tends to both reduce number density (where there is rotation) and increase it (where there is strain), and eventually, the number density tends to a distribution that is mostly independent of initial condition.

An external file that holds a picture, illustration, etc.
Object name is nihms864032f6.jpg

Fig. 6

Insensitivity to the initial condition. For panels (a) and (b) droplets were placed randomly on each 1 m × 1 m square. For (c) and (d) droplets were arranged on a regular grid on each 50 cm × 50 cm square.

c. Pulse to pulse fluctuation and averaging

The SNR1 results in §3a were based on using the first term in Equation (23), which sums the powers reflected by individual droplets. It was argued that for a statistically stationary target, this should equal the average power from many pulses. In the present case, the droplet configuration is not spatially homogeneous and is descending through a fixed beam. Hence, the question arises whether the powers returned from a sequence of pulses can be considered to be statistically stationary in a certain interval, and if so, how many pulses is sufficient to recover the SNR1 values presented.

To obtain complex voltage returns from a sequence of pulses one needs to evolve the wake in time, however, the method that was described in §2 gives a trail of droplets at a single instant of time, t. To evolve this configuration to time t + Δt, the configuration at t is translated horizontally by Δx = −U_appΔt, i.e., the droplet trail is assumed to be invariant in a reference frame moving to the left with the airplane. This procedure does not correspond exactly to reality, but captures both the rotation of droplets around the vortices, and their vertical descent with time at a fixed location. The received complex voltage is evaluated using (18) at a sequence of times separated by the pulse repetition period, keeping the resolution volume centered at (x,y,z) = (6 nm,−50 m,−230 m). This location corresponds to the upper SNR1 peak of the crescent wrench in Figure 4. The value of the pulse repetition frequency (PRF) was chosen to be at or close to the highest value available for each radar.

Figure 7 shows the result for a time period during which a cluster of droplets enters and leaves the beam. The power in individual pulse returns is shown in gray. The total number of active pulses changes from radar to radar because their beam widths and PRFs are different. In particular, the period of activity was found to equal the time it would take the vortex pair to descend through roughly one-third of the vertical projection of the half-power beam width. The average of pulse powers is shown in green over an averaging segment whose length is 512 pulses. The red curve shows the value of SNR1. Our assumption was that SNR1 should equal the green level. This is seen to be true to a good degree. The fluctuations are due to statistical error and were found to decrease with increasing the averaging interval. It is worth remembering here that the radar has access to only the individual pulse returns and their average (for example the green values); only the simulation has access to the red curve (SNR1).

An external file that holds a picture, illustration, etc.
Object name is nihms864032f7.jpg

Fig. 7

Gray line: Instantaneous power P_r(t) received from the same range gate due to a sequence of transmitted pulses. The range cell is at (x,y,z) = (6 nm,−50 m,−230 m). Green: received powers averaged over segments 512 pulses long. Red: the first term in (23). Panels (a)–(e) are for the same cases as in Figure 4b–f.

d. Convergence of pulse statistics

Recall that each computational droplet at a single location represents M_true actual droplets located at different positions; in fact M_true≈ 100 in the calculations presented. It was claimed (§2g) that this should not affect pulse statistics, provided the number of computational droplets is sufficiently large. To verify this, the number of computational droplets was increased by four. Random placement of droplets was employed in the injected squares. Four realizations of the droplet trail were generated using different random number seeds for the initial size distribution and droplet placement. The four realizations were then merged into one trail for the radar reflectivity calculation. Figure 8 shows that the probability density p(|V|) of the modulus |V| of complex voltage is unchanged by the resolution refinement. The probability densities are very well fit by the Rayleigh distribution (Beckmann 1962)

An external file that holds a picture, illustration, etc.
Object name is nihms864032f8.jpg

Fig. 8

Convergence of pulse statistics when the number of computational droplets is increased by four. The same case as Figure 7d is used (MIRA-35 radar). Note: The solid and dashed lines are nearly coincident.

p (| V |) = \frac{| V |}{σ_{R}^{2}} exp (- {| V |}^{2} / 2 σ_{R}^{2}),

(40)

having the same mean as the data. The Rayleigh distribution results when the scattering amplitude is the same for all droplets and the phases uniformly distributed. The case selected is the same as that shown in Figure 7d (apart from the random placement of droplets in the injected squares). Pulses in the interval of stationarity were chosen, namely pulse number ∈ [−8000,8000].

e. A non-IFR condition

It would be valuable to have the capability to detect wakes in non-IFR conditions. Furthermore, in a flight test study of the feasibility of the present proposal, it would be too costly to wait until IFR conditions occur before a test can be conducted. For this reason it is of interest to know what reflectivity is obtained at less humid and less cold conditions. We chose RH = 60% and T = 20 C.

With the previous choice of n_square = 15 and n_x = 120 as injection parameters, it was found that a high rate of evaporation resulted in a small number of computational droplets remaining near x = 6 nm. This increased statistical error. To reduce sampling error, an ensemble of ten trails were computed with different random number seeds for the droplet size sample. The ensemble was then combined into a single trail for the reflectivity analysis. As a result, the total number of computational droplets is so large that each one presents only 9.8 true droplets in the reflectivity analysis.

Figure 9 shows that only the high-frequency radars, MIRA-35 and W-SACR, give positive values of SNR1 (dB) in the vicinity of the vortices, and even these values are marginal. To increase SNR1, the number of nozzles could be increased; for instance four nozzles on each side of the aircraft would increase SNR1 by 6 dB.

An external file that holds a picture, illustration, etc.
Object name is nihms864032f9.jpg

Fig. 9

SNR1 for a non-IFR condition (RH = 60%, T = 20 C). A range elevation scan of the wake cross-section 6 nm behind the aircraft is shown. The white circles in panel (e) are points for which a spectral analysis is presented in Figure 10.

There is a powerful method that enables detection even when SNR1 (dB) < 0. It comes at the cost of increased dwell and processing time. We learnt about the method from notes on the sensitivity of the MIRA-35 radar given to us by Matthias Bauer-Pfundstein (Metek). It is also briefly described in Görsdorf et al. (2015, p 680). The idea is that in a discrete Fourier transform, the noise is spread equally to all the frequency bins, whereas the spectrum of the signal is confined to only a few of the bins. (The latter is true provided the probability distribution of droplet velocities in the resolution cell is narrow compared to 2U_max for the radar. For MIRA- 35, for example, at PRF = 10 kHz we have 2U_max = 42 m s⁻^{and so this is unlikely to be an issue.) Hence, an FFT effectively reduces the noise by a factor of}N_FFT.

To investigate this technique, complex white noise with a mean power equal to P_noise for the radar was added to complex voltages of pulse returns. Illustrative results are shown in Figure 10. Panels (a) and (b) are for a range cell centered at the left white dot in Figure 9e where SNR1 = −2.7 dB. Panels (c) and (d) are for the right white dot where SNR1 is even lower, namely, −7.2 dB. Consider panels (a) and (b). An averaging of pulse power returns by the radar would give values (the green line) only slightly above the noise, not enough for a positive detection. Averaging the doppler spectra from 10 segments gives panel (b) with a peak 40 dB above the noise. In the present example, this would require a dwell time of 0.5 s for each elevation angle. For the second location where SNR1 is weaker, the doppler spectrum has a peak that is about 25 dB above the noise (using the same dwell time).

An external file that holds a picture, illustration, etc.
Object name is nihms864032f10.jpg

Fig. 10

Detection at low SNR1 using spectral processing. Panels (a) and (b) are for the the resolution cell centered on (x,y,z) = (6 nm,−55 m,−260 m) which is shown as the white circle to the left in Figure 9e. Panels (c) and (d) are for (x,y,z) = (6 nm,−15 m,−255 m) which is shown as the white circle to the right in Figure 9e. The radar is MIRA-35 with τ = 0.2 μs and PRF = 10 kHz. Non-IFR condition (RH = 60%, T = 20 C).

f. Power-weighted average radial velocity

It has been stated (Doviak and Zrnić 1984, §5.2) that the first moment of the doppler spectrum is the radial velocity of droplets in the resolution volume, weighted by their individual scattered powers. This is reasonable although we have neither encountered nor attempted a proof of it starting from Equation (2). Figure 11 displays the power-weighted average radial velocity corresponding to the cases previously shown in Figure 4. Only points where SNR (dB) > 0 are shown. The actual radial velocity (with respect to the radar) of the gas is shown in panel (f). The radar data appears as a filtered version of the actual velocity and, due to particle centrifugation, is unable to detect the maximum value of 19.8 m s⁻^{in the vortex core. Nevertheless, the radars give a good representation of the gas velocity where particles are present. To estimate vortex circulation, the observed velocity is multiplied by 2}πr where r is the distance from the vortex center and can be determined from the location of zero radial velocity.

An external file that holds a picture, illustration, etc.
Object name is nihms864032f11.jpg

Fig. 11

Power-weighted radial velocity corresponding to Figure 4. Only points where SNR1 > 1 are colored. The actual radial velocity of the gas is shown in panel (a).

g. Effect of vortex core growth

The vortex model presented in §2f assumed that the vortex core radius does not grow with downstream distance. In reality, the core radius grows and the peak velocity diminishes. Experiments on aircraft wakes (see Govindaraju and Saffman 1971) have shown that the core radius, r₁, defined to be the radius where the tangential velocity peaks, diffuses as follows with distance x behind the plane:

r_{1} = b_{1} {[\frac{(x - x_{0}) Γ}{U_{app}}]}^{1 / 2},

(41)

where b₁ is a constant that depends on the vortex Reynolds number, Γ/ν. In the present case Γ/ν = 3.5 × 10^{, or which Table I in}Govindaraju and Saffman (1971) gives b₁ = 1.3 × 10⁻^{(using the corrected value listed in their table). The value of the parameter}x₀ (called the virtual origin) can be obtained from the core radius at x = 0. Plotting the profile we find that r₁ = 0.8 m for the present case which gives x₀ = −558.8 m.

To include core growth in the vortex model of §2f, the variable η is redefined as

η \equiv \frac{r}{g (x) b_{0}},

(42)

where g is the growth function

g (x) = \frac{r_{1} (x)}{r_{1} (0)} = {(\frac{x - x_{0}}{- x_{0}})}^{1 / 2} .

(43)

Figure 12 displays the decay of peak tangential velocity from x = 0 to x = 6 nm.

An external file that holds a picture, illustration, etc.
Object name is nihms864032f12.jpg

Fig. 12

Decay of tangential velocity due to core diffusion.

Figures 13a and b show the droplet configuration and SNR1, respectively, for the MIRA-35 radar when core growth is included; they should be compared with Figures 4a and e. This comparison shows that slightly less centrifugation takes places with core diffusion which makes the region of reflectivity slightly smaller.

An external file that holds a picture, illustration, etc.
Object name is nihms864032f13.jpg

Fig. 13

Effect of vortex core diffusion. (a) Droplet configuration. (b) SNR1 in the x = 6 nm cross-section. IFR ambient conditions (RH = 92.7%, T = 15.2 C). Nozzle 1.

h. Comparison of four nozzles

Finally, Table 6 compares SNR1 reflectivity for the four nozzles listed in Table 2. The range cell for all cases is centered at (x,y,z) = (6nm,−50m,−230m), which corresponds to the top of the crescent wrench in Figure 4. IFR ambient conditions have been assumed. It is observed that the quantity ${\dot{ζ}}_{nozzle}^{<}$ (defined in Equation 9), which depends only on the droplet size distribution produced by a nozzle, provides an excellent indicator of the relative performance of different nozzles.

TABLE 6

A comparison of SNR1 reflectivity obtained from a single range cell by using the four different nozzles listed in Table 2. The cell center is at (x,y,z) = (6 nm,−50 m,−230 m). The signed values (±) indicate values relative to nozzle 1. MIRA-35, τ = 0.2 μs, θ_b = 0.52°.

Nozzle no.	gpm	SNR1 (dB)	˙ζ< (dB m^s⁻^{) nozzle}	ζ_x (dB m^{) near vortex at}x = 6 nm
1	3.70	19.15	−157.4	−178.8
2	6.49	+1.03	+1.0	+0.6
3	3.00	−1.04	−1.2	−1.4
4	3.06	−0.5	−0.3	−1.0

To provide further insight we compute the quantity ζ_x which is defined to be ζ (see Equation 8) per unit axial length of the wake. It is calculated as a diagnostic of the droplet trajectory and size evolution and is shown in Figure 14. The solid lines give the total value (over an entire cross-section) and diagnose total evaporative loss. The dashed lines give the contribution from droplets in a neighborhood (defined in the caption) of the vortex: these curves diagnose both evaporative loss and loss by sedimentation. The ordering of ζ_x values (pertaining to the neighborhood of the vortex) at x = 6 nm, which are also listed in the last column of Table 6, matches the ordering of SNR1 for the different nozzles.

An external file that holds a picture, illustration, etc.
Object name is nihms864032f14.jpg

Fig. 14

(a) ζ_x for a droplet trail on one side of the aircraft. “Near the vortex” curves (dashed) were obtained by considering only those droplets that obey |y| < 60 m and |z−z_vort| < 35 m, z_vort being the height of the vortex center. IFR ambient humidity and temperature were assumed. (b) Droplet size distributions produced by the four nozzles. The solid lines show the exact log-normal distribution, while the symbols show the distribution for each sample of 27,000 droplets injected into the wake.

All statements of comparison in the following are relative to nozzle 1 and make reference to Figure 14a. If the increased SNR1 reflectivity of nozzle 2 were due to increased volume alone, we would get a 2.4 dB increase in reflectivity. The actual increase is 1.03 dB. To understand this, we first observe that nozzle 2 (solid red line) initially has a 0.7 dB higher value of ζ_x, less than the 2.4 dB increase in its volume flow-rate. This is because nozzle 2 produces more small droplets. By x = 6 nm the 0.7 dB increase has been reduced to 0.55 dB because the smaller droplets of nozzle 2 evaporate faster. The fact that ζ_x in the vicinity of the vortex is 0.6 dB higher must arise from the fact that the smaller droplets of nozzle 2 have sedimented less.

Consider nozzle 3 (green curves). If its smaller reflectivity (relative to nozzle 1) were due to decrease in volume, then we would expect a −0.91 drop in SNR1 which is close to what is obtained. This is understandable given that its initial ζ_x is very nearly the same as for nozzle 1. This is surprising given that nozzle 3 has many more smaller droplets. However, close inspection of its size distribution (green line in Figure 14b) shows that it also has more droplets that are very large (specifically a > 370 μm). This fact also explains the more rapid loss of ζ_x by sedimentation (dashed green curve in Figure 14a) and less rapid loss by evaporation (solid green curve). Overall, these two effects balance and the final effect that remains is that due to volume decrease.

Despite its smaller flow-rate, nozzle 4 has a higher initial value of ζ_x; see the solid blue curve. This is because it produces more large droplets. Unfortunately, they rapidly fall out of the wake (dashed blue curve).

a. Signal-to-noise ratio for the IFR case

Figure 5 shows that with its lowest available pulse width of 0.05 μs, the W-SACR radar is able to resolve some of the spiral structure of the droplet pattern at the expense of some loss in SNR1.

Fig. 5

Simulated SNR1 for W-SACR with a pulse width of 0.05 μs.

p (| V |) = \frac{| V |}{σ_{R}^{2}} exp (- {| V |}^{2} / 2 σ_{R}^{2}),

(40)

e. A non-IFR condition

Fig. 9

Fig. 10

f. Power-weighted average radial velocity

Fig. 11

Power-weighted radial velocity corresponding to Figure 4. Only points where SNR1 > 1 are colored. The actual radial velocity of the gas is shown in panel (a).

g. Effect of vortex core growth

r_{1} = b_{1} {[\frac{(x - x_{0}) Γ}{U_{app}}]}^{1 / 2},

(41)

To include core growth in the vortex model of §2f, the variable η is redefined as

η \equiv \frac{r}{g (x) b_{0}},

(42)

where g is the growth function

g (x) = \frac{r_{1} (x)}{r_{1} (0)} = {(\frac{x - x_{0}}{- x_{0}})}^{1 / 2} .

(43)

Figure 12 displays the decay of peak tangential velocity from x = 0 to x = 6 nm.

Fig. 12

Decay of tangential velocity due to core diffusion.

Fig. 13

Effect of vortex core diffusion. (a) Droplet configuration. (b) SNR1 in the x = 6 nm cross-section. IFR ambient conditions (RH = 92.7%, T = 15.2 C). Nozzle 1.

h. Comparison of four nozzles

TABLE 6

Nozzle no.	gpm	SNR1 (dB)	˙ζ< (dB m^s⁻^{) nozzle}	ζ_x (dB m^{) near vortex at}x = 6 nm
1	3.70	19.15	−157.4	−178.8
2	6.49	+1.03	+1.0	+0.6
3	3.00	−1.04	−1.2	−1.4
4	3.06	−0.5	−0.3	−1.0

Fig. 14

4. Concluding Remarks

It was proposed that spraying a small amount of water into the vortex wake of a heavy aircraft during landing can make the wake visible to existing weather/cloud radars and thereby aid air traffic controllers in selecting appropriate aircraft separations. This approach could also be used for wake vortex studies of aircraft.

Simulations of the radar reflectivity of the spray trail were performed for existing weather/cloud radars. For ambient humidity at the lower end of values typical for IFR conditions, the results showed that that good signal-to-noise (SNR) ratios (averaged over many pulses) are obtained at distances behind the aircraft of up to 6 nm, the largest that would be contemplated given existing wake separations used in air traffic control. For the case most studied here, the amount of water spray was 3 gallons per nautical mile of wake that needs to be detected. A currently available nozzle used for agricultural spraying can be used. A doubling of volume by doubling the number of nozzles gives a proportional increase in SNR. For a case of average humidity, evaporation for severe and pulse-averaged, SNR values dropped below unity. However, since the pulse returns of the wake remained statistically stationary for 1 to 6 secs (depending on the radar), it was shown that the signal-to-noise ratio can be increased to detectable levels by spectral (doppler) processing and averaging doppler spectra for consecutive time segments. This would require greater dwell time for each direction the radar is pointed at. Ultimately, selecting the dwell time for a given situation will be a trade-off between quickly completing a scan of a wake cross-section and increasing SNR.

1. Suggested future work

As an airplane nears the touch-down point, flaps are deflected at increasing angles. The presence of flap vortices should be included in future analysis.
The present work has ignored space-time fluctuations of the air velocity field. They will arise from the direct effect of atmospheric turbulence and from vortex core waviness induced by atmospheric turbulence, and further amplified by vortex core instabilities. Velocity fluctuations will disperse the spray trail and if this happens on the scale of the pulse width or beam width, then reflectivity will be reduced. This effect should be studied in future work.

2. Application notes

For the purposes of simulation we generated a spray trail that was 7 nm long. In practice, to reduce the volume of water, spray would be released only at axial locations where a detection would be performed. For each detection location, the length of the trail would need to be a few beam widths long and the release location would have to account for any head/tail wind. A trail that is three beam widths long would require only 0.084 gallons. This value assumes that θ_b = 1°, range = 1 km, and a flow-rate of 3 gallons nm⁻^{(counting both sides of the airplane). Hence, there is considerable room for increasing water volume, and therefore signal-to-noise ratio. The main difficulty is that for the nozzles presently considered, more than one would be required. A better solution might be to design a spray head containing several nozzles.}
Given that spectral processing is required for detection in conditions of average humidity, it is likely that processing decisions will have to be based on humidity or on the quality of incoming returns. If the humidity is high and the quality of returns high, then the mean velocity can be obtained from a pulse-pair estimate. If the humidity is low, then spectral processing can be turned on.
Some of the requirements of the present application are similar to those for radar imaging of tornados (French et al. 2014). This includes a smaller detection volume and the need to complete a scan faster than the vortex evolution time. Therefore the technology developed for that application could be useful here.
In IFR conditions, natural precipitation (fog, mist, drizzle, or heavy rain) will be present between the radar and the wake and lead to absorption. However, at ranges of ≈ 1 nm envisioned for the present application, this is small.
If spraying is to be employed in very cold conditions (Denver, Colorado comes to mind), freezing of water must obviously be prevented in the water storage and delivery system.
Dual polarization. Droplets moving relative to the air become oblate due to a higher air pressure at the front stagnation point and low pressure at 90° from the front stagnation point. For falling rain droplets, this results in greater reflected power from incident waves that are horizontally versus vertically polarized (Doviak and Zrnić 1984, §8.5.3). Most weather radars employ dual polarization to obtain more information about rainfall rate. Since, in the present case, droplets revolving around the vortices are small and their velocity relative to the air is also small, we expect that droplets will remain very nearly spherical. Therefore, it is not expected that dual polarization would provide additional information about the flow. However Keränen and Chandrasekhar (2014) have suggested that dual polarization could be used for enhancing SNR. This works by exploiting coherence between signals in the horizontal and vertical channels.
Since the maximum range pertinent to the present application is much lower than for cloud and precipitation detection, the pulse repetition frequency could be increased (the maximum duty cycle of the klystron or magnetron permitting) in order to reduce the dwell time for spectral averaging.
One obvious modification of existing cloud/precipitation radar software for the present application would be a reduction in the spacing of range gates from their current values, for example 25 m which is employed in Ka-SACR and W-SACR (Kollias et al. 2014).

1. Suggested future work

As an airplane nears the touch-down point, flaps are deflected at increasing angles. The presence of flap vortices should be included in future analysis.
The present work has ignored space-time fluctuations of the air velocity field. They will arise from the direct effect of atmospheric turbulence and from vortex core waviness induced by atmospheric turbulence, and further amplified by vortex core instabilities. Velocity fluctuations will disperse the spray trail and if this happens on the scale of the pulse width or beam width, then reflectivity will be reduced. This effect should be studied in future work.

2. Application notes

For the purposes of simulation we generated a spray trail that was 7 nm long. In practice, to reduce the volume of water, spray would be released only at axial locations where a detection would be performed. For each detection location, the length of the trail would need to be a few beam widths long and the release location would have to account for any head/tail wind. A trail that is three beam widths long would require only 0.084 gallons. This value assumes that θ_b = 1°, range = 1 km, and a flow-rate of 3 gallons nm⁻^{(counting both sides of the airplane). Hence, there is considerable room for increasing water volume, and therefore signal-to-noise ratio. The main difficulty is that for the nozzles presently considered, more than one would be required. A better solution might be to design a spray head containing several nozzles.}
Given that spectral processing is required for detection in conditions of average humidity, it is likely that processing decisions will have to be based on humidity or on the quality of incoming returns. If the humidity is high and the quality of returns high, then the mean velocity can be obtained from a pulse-pair estimate. If the humidity is low, then spectral processing can be turned on.
Some of the requirements of the present application are similar to those for radar imaging of tornados (French et al. 2014). This includes a smaller detection volume and the need to complete a scan faster than the vortex evolution time. Therefore the technology developed for that application could be useful here.
In IFR conditions, natural precipitation (fog, mist, drizzle, or heavy rain) will be present between the radar and the wake and lead to absorption. However, at ranges of ≈ 1 nm envisioned for the present application, this is small.
If spraying is to be employed in very cold conditions (Denver, Colorado comes to mind), freezing of water must obviously be prevented in the water storage and delivery system.
Dual polarization. Droplets moving relative to the air become oblate due to a higher air pressure at the front stagnation point and low pressure at 90° from the front stagnation point. For falling rain droplets, this results in greater reflected power from incident waves that are horizontally versus vertically polarized (Doviak and Zrnić 1984, §8.5.3). Most weather radars employ dual polarization to obtain more information about rainfall rate. Since, in the present case, droplets revolving around the vortices are small and their velocity relative to the air is also small, we expect that droplets will remain very nearly spherical. Therefore, it is not expected that dual polarization would provide additional information about the flow. However Keränen and Chandrasekhar (2014) have suggested that dual polarization could be used for enhancing SNR. This works by exploiting coherence between signals in the horizontal and vertical channels.
Since the maximum range pertinent to the present application is much lower than for cloud and precipitation detection, the pulse repetition frequency could be increased (the maximum duty cycle of the klystron or magnetron permitting) in order to reduce the dwell time for spectral averaging.
One obvious modification of existing cloud/precipitation radar software for the present application would be a reduction in the spacing of range gates from their current values, for example 25 m which is employed in Ka-SACR and W-SACR (Kollias et al. 2014).

Acknowledgments

Acknowledgments. I am grateful to several individuals for sending me information when requested and sometimes more information than requested. Matthias Bauer-Pfundstein (Metek) sent me information on the MIRA-35 radar including a detailed dBZ sensitivity calculation, information on signal processing, and detailed information on loss determination. It was from his dBZ sensitivity document that I learnt that a discrete Fourier transform enhances signal-to-noise ratio in the spectral bin of interest. John Cho (MIT Lincoln Laboratories) sent me information on the TDWR (Terminal Doppler Weather Radar). Brad Fritz (U.S. Dept. of Agriculture) sent me his Excel program giving parameters of the droplet size distribution produced by different aerial spray nozzles at different free-stream air speeds. He also experimentally measured for me the flowrate of the Davidon-Triset nozzle. Calvin Kroes (CP Products) sent me information on nozzle orifice areas, flow rates, and exit velocities for various CP nozzles. Keith Vickers (Enterprise Electronics Corporation) sent me basic information on the DWSR series of radars. I am grateful to Alan Wray (NASA Ames) for useful discussions. I am grateful to Jasim Ahmed and Alan Wray for performing the internal review.

A1. Drag coefficient

With the Reynolds number defined as Re ≡ 2a|u_rel|/ν_a (where ν_a is the kinematic viscosity of air), the coefficient of drag, C_D, is obtained by numerically inverting the relation (PK, Equation 10–145):

Y = \sum_{m = 0}^{6} B_{m} X^{m},

(A1)

where X = ln(C_DRe^{) and Re = exp(}Y). This formula is originally from Beard (1976) and is based on the drag coefficient of a solid sphere. The validity of this rests on two assumptions. The first is that the droplet does not distort significantly from being spherical. The equilibrium aspect ratio of a falling raindrop is given by PK (Equation 10–108):

\frac{b}{a} = \frac{1 - 0.11 We}{1 + 0.11 We} .

(A2)

Here $We \equiv 2 a ρ_{a} u_{rel}^{2} / γ_{w / a}$ is the Weber number, where γ_w_/_a is the surface tension of water in air. For nozzle 1 and the IFR case, the smallest value of b/a was 0.9 which occurred at early times for a droplet which quickly fell below the vortex. We conclude that droplet deformation is negligible particularly for those droplets that remain with the vortex pair. The second assumption is that the ratio of dynamic viscosities, η_a/η_w≈ 1.8 × 10⁻^{is small. In the creeping flow limit, the Hadamard-Rybczynski formula (see}Beard 1976) for the drag of a water sphere divided by the drag of a solid sphere is

\frac{F_{D}}{F_{D s}} = 1 - η_{a} / 3 η_{w} = 0.995.

(A3)

Numerical solutions (PK, p. 388) indicate that for Re < 300, the drag coefficient of a water sphere differs by less than ∼ 1% from that of a solid sphere.

Formula (A1) is valid for Re < 500. Only for nozzle 4 was this condition slightly exceeded for a few droplets. For Re ≤ 1.5, the explicit formula for solid spheres (White 1974, eq. 3–265) was found to agree well with (A1) and was used instead. We have also implemented but not used the Schiller and Naumann explicit drag formula for a solid sphere (e.g., Apte et al. 2003)

C_{D} = \frac{24}{Re} ((1 + 0.15 {Re}^{0.687})

(A4)

which is said to be accurate to within 5% for Re < 800. Figure 15 compares the three formulas for C_D up to the maximum value of Re = 800 we allow in the code. It suggests that in the future it would be as accurate to use the explicit Schiller and Naumann formula, which is cheaper to compute.

An external file that holds a picture, illustration, etc.
Object name is nihms864032f15.jpg

Fig. 15

Comparison of three formulas for the drag coefficient of a solid sphere.

A2. Evaporation model

The evaporation model uses appropriate formulas from PK. These formulas are collected here to document the choices we have made and also because they are scattered throughout the book. Gas kinetic effects and the Kelvin curvature effect are neglected since we remove droplets when their radii fall below 20 μm. Throughout, T_C denotes centigrade temperature:

T_C = T − 273.15.

(A5)

Note that there are 100 Pa in a mb. Subscripts ‘a’ and ‘w’ denote air and water, respectively.

a. Properties of air, liquid water, and water vapor

Given the universal gas constant, ℛ = 8.3144 J K^mol⁻^{, and the molar mass of dry air,}M_a = 28.9644 × 10⁻^{kg mol}⁻^{, the gas constant for air is}

R_a = ℛ/M_a.

(A6)

The density of air is calculated from the ideal gas law:

ρ_a = p/R_aT.

(A7)

The density of liquid water at p = 1 atm is given by (PK 3–13):

10^{3} ρ_{w} = \frac{\sum_{m = 0}^{5} A_{m} {T_{C}}^{m}}{1 + B T} g {cm}^{- 3}, 0 \leq T_{C} \leq 100,

(A8)

with A₀ = 999.8396,A₁ = 18.224944,A₂ = −7.922210 × 10⁻^,A₃ = −55.44846 × 10⁻^,A₄ = 149.7562 × 10⁻^,A₅ = −393.2952 × 10⁻^,B = 18.159725 × 10⁻^{. The thermal diffusivity of air is}

κ_a = k_a/ρ_aC_p,

(A9)

where the conductivity is given by (PK 13–18a):

k_a = (5.69 + 0.017T_C) × 10^{cal cm}^s^K^.

(A10)

Note that there are 4.184 J per cal. The heat capacity of air is:

C_p = 1006.1 J kg^K^.

(A11)

The diffusivity, D_v, of water vapor is calculated using Equation (PK 13–3):

D_{v} = 0.211 {(\frac{T}{T_{0}})}^{1.94} (\frac{p_{0}}{p}) {cm}^{2} s^{- 1},

(A12)

with T₀ = 273.15 K and p₀ = 1013.25 mb. The dynamic viscosity of air is (PK 10–141)

η_{a} = (1.718 + 0.0049 T_{C}) \times 10^{- 4} poise, T_{C} \geq 0.

(A13)

The Schmidt number of vapor is defined (PK, p. 538) as

Sc_v ≡ ν_a/D_v,

(A14)

where ν_a = η_a.ρ_a is the kinematic viscosity. The Schmidt number for heat is

Sc_h ≡ ν_a/κ_a.

(A15)

b. Evolution of droplet radius

The evolution of droplet radius a(t) is given by:

a \frac{d a}{d t} = {(a \frac{d a}{d t})}_{0} f_{v},

(A16)

where f_v, called the ventilation coefficient, represents the enhancement of evaporation rate due to advection of air past the droplet, and ()₀ represents a quantity in the absence of advection.

Let Re = 2aU_rel/ν_a denote the Reynolds number based on droplet diameter and drop speed U_rel relative to the air. Defining $F \equiv {Sc}_{v}^{1 / 3} {Re}^{1 / 2},$ the ventilation coefficient is given by (PK 13–60) and (PK 13–61):

f_{v} = {\begin{cases} 1.00 + 0.108 F^{2}, & F < 1.4; \\ 0.78 + 0.308 F, & 1.4 \leq F < 51.4. \end{cases}

(A17)

The first factor on the RHS of (A16), which represents evaporation in the absence of advection, is given by

{(a \frac{d a}{d t})}_{0} = \frac{D_{v} M_{w}}{ℛ ρ_{w}} (\frac{e_{\infty}}{T_{\infty}} - \frac{e_{sat} (T_{a})}{T_{a}}),

(A18)

where D_v is the vapor diffusivity (which we evaluated at ambient conditions using equation A12), M_w = 28.97 gm mol⁻^{is the molecular mass of water,}e_∞ is the vapor pressure in the ambient, T_∞ is the ambient temperature, and e_sat(T_a) is the saturation vapor pressure evaluated at the surface temperature T_a of the droplet. For e_sat(T) we use the expression (Sonntag 1994, eq. 7)

e_{sat} (T) = 100 exp (\sum_{m = 1}^{4} a_{m} T^{m - 2} + a_{5} ln T) Pa, 173.15 \leq T \leq 373.15,

(A19)

with a₁ = −6.0969385 × 10^,a₂ = 1.6635794 × 10^,a₃ = −2.711193 × 10⁻^,a₄ = 1.673952 × 10⁻^,a₅ = 2.433502.

c. Temperature at the droplet surface

The internal energy of the water droplet is:

q = mC_wT_a,

(A20)

where m is its mass, C_w = 4.187 × 10^{J kg}⁻^K⁻^{is the heat capacity of water, and}T_a is its temperature, which we have taken to be uniform and equal to the value at the surface. The internal energy of the drop increases due to diffusion of heat at its surface and release of latent heat (PK 13–65):

\frac{d q}{d t} = 4 π f_{h} a k_{a} (T_{\infty} - T_{a}) + L_{e} \frac{d m}{d t},

(A21)

where L_e is the latent enthalpy of evaporation of pure water evaluated at the surface temperature of the drop. The dependence of L_e on centigrade temperature T_C is:

L_e = (2500.8 − 2.36T_C + 0.0016T_C^{− 0.00006}T_C^{) × 10}^{J kg}^.

(A22)

Substituting (A20) into (A21) gives

\frac{d T_{a}}{d t} \frac{3 f_{h} k_{a} \infty}{C_{w} a^{2} ρ_{w}} (T_{\infty} - T_{a}) + 3 (L_{e} / C_{w} - T_{a}) \frac{1}{a} \frac{d a}{d t} .

(A23)

Here f_h is the ventilation coefficient for heat: it is given by the same expression as (A17) except with $F = {Sc}_{h}^{1 / 3} {Re}^{1 / 2},$ , where Sc_h is the Schmidt number for heat.

Karim Shariff, NASA Ames, Moffett Field, California, 94035;

^{Corresponding author address: NASA Ames, Moffett Field, California, 94035,}vog.asan@ffirahS.miraK

Abstract

Aircraft trailing vortices pose a danger to following aircraft during take-off and landing. This necessitates spacing rules, based on aircraft type, to be enforced during approach in IFR (Instrument Flight Regulations) conditions; this can limit airport capacity. To help choose aircraft spacing based on the actual location and strength of the wake, it is proposed that wake vortices can be detected using conventional precipitation and cloud radars. This is enabled by spraying a small quantity water into the wake from near the wing. The vortex strength is revealed by the doppler velocity of the droplets. In the present work, droplet size distributions produced by nozzles used for aerial spraying are considered. Droplet trajectory and evaporation in the flow-field is numerically calculated for a heavy aircraft, followed by an evaluation of radar reflectivity at 6 nautical miles behind the aircraft. Small droplets evaporate away while larger droplets fall out of the wake. In the humid conditions that typically prevail during IFR, a sufficient number of droplets remain in the wake and give good signal-to-noise ratios (SNR). For conditions of average humidity, higher frequency radars combined with spectral processing gives good SNR.

Abstract

Footnotes

^{A former colleague, Dr. Vernon Rossow, has long suggested that with GPS and fly-by-wire, the simplest approach to wake avoidance would be for a following aircraft to remain at or above the path of the leader.}

Footnotes