## Abstract

Simulation of radar beam propagation is an important component of numerous radar applications in meteorology, including height assignment, quality control, and especially the so-called radar forward operator. Although beam propagation in the atmosphere depends on the refractive index and its vertical variation, which themselves depend on the actual state of the atmosphere, the most common method is to apply the 4/3 earth radius model, based on climatological standard conditions. Serious deviations from the climatological value can occur under so-called ducting conditions, where radar beams at low elevations can be trapped or propagate in a waveguide-like fashion, such that this model is unsuitable in this case. To account for the actual atmospheric conditions, sophisticated methods have been developed in literature. However, concerning the practical implementation of these methods, it was determined that the description in the literature is not always complete with respect to possible pitfalls for practical implementations.

In this paper, a revised version of an existing method (one example for the above-mentioned “pitfall” statement) is introduced that exploits Snell’s law for spherically stratified media. From Snell’s law, the correct sign of the local elevation is a priori ambiguous, and the revised method explicitly applies (i) a total reflection criterion and (ii) another ad hoc criterion to solve the problem.

Additionally, a new method, based on an ordinary differential equation with respect to range, is proposed in this paper that has no ambiguity.

Sensitivity experiments are conducted to investigate the properties of these three methods. The results show that both the revised and new methods are robust under nonstandard conditions. But considering the need to catch an elevation sign ambiguity in the revised method (which cannot be excluded to fail in rare instances), the new method is regarded as more robust and unproblematic, for example, for applications in radar forward operators.

## 1. Introduction

The microwave pulses emitted by any radar system do not propagate through the atmosphere along straight lines but are usually more or less curved toward the earth’s surface because of atmospheric refraction. Refraction is caused by the (more or less gradual, or “bumpy”) decrease of the real part of the air refractive index *n* from values >1 near the ground to 1 in outer space. Here, especially the gradient component of *n* perpendicular to the travel direction of the radiation is important.

There are a number of applications that require a more or less precise knowledge about this so-called radar beam propagation in the atmosphere, and these are all applications where the beam height above the earth’s surface as a function of antenna elevation *ϵ*_{0} and distance from the radar is crucial. The meteorological applications range from the correct height assignment of radar echoes as function of range to the monitoring of “anomalous” propagation conditions (with their associated phenomena-like abnormal ground clutter, radar “blind” zones, and dispersion lens effects) in radar data quality control.

Another application is in the so-called radar forward operators for “synthetic” radars, which simulate the measurement process and outputs of a radar from prognostic model outputs and allow for comparisons in terms of direct radar data–like reflectivity and Doppler velocity, for the purposes of numerical weather prediction (NWP) model validation and data assimilation. One part of the radar simulation process is the computation of beam propagation within the atmosphere simulated by an NWP model in an appropriate way. Low elevations are often vulnerable to anomalous propagation and orographic beam blockage, which can seriously affect the radar’s ability to detect and quantify precipitation at ground level. Important issues here are to minimize the influences of these effects in the observation–simulation comparison.

Much of the fundamental research on radar wave propagation in the atmosphere has been stimulated by military demands during World War II. Comprehensive summaries and textbooks have been published on the subject shortly after the war, for example, Kerr (1951), among others. Most of the theoretical background known at that time is still state of the art. Based on these established theories, the problem of radar wave propagation in radar meteorology is usually reduced to the tracing of single rays according to classical geometrical optics (ray tracing).

Practical approximations and algorithms have been derived, and as of today, there are several approaches in the meteorological literature. The most popular one is the approximate ⅘ earth radius model by Schelleng et al. (1933), which is widely used in radar meteorology because of its simple conceptual formulation and relative ease of implementation. It has, however, limitations in nonnormal propagation conditions because it only takes average climatological atmospheric conditions into account. Already Hartree et al. (1946) derived a differential equation for ray heights as function of surface distance based on actual atmospheric conditions and solved an approximated variant of it. Freehafer (1951) presented the theoretical derivation of an integral conserved quantity of this equation (Snell’s law in a spherically stratified medium) but did not make the connection between the two. It has to be stated that the above-mentioned variant of Snell’s law was known to astronomers and meteorologists much earlier, see, for example, Perntner and Exner (1922). Today, almost all methods in the literature that take into account the actual atmospheric conditions are iterative methods based on either the discretized variants of Snell’s law (e.g., Caumont et al. 2006; Chen et al. 2009) or generalizations of the ⅘ earth radius model (e.g., Gao et al. 2006). However, concerning the practical implementation of these methods, we find the description in the literature sometimes does not complete and satisfy with respect to possible pitfalls (singularities, ambiguities, etc.).

In the present paper, we discuss the pitfalls of the method of Caumont et al. (2006; in this case, it is how to determine the sign of the ray’s local elevations), give a solution to the problem, and revise the method accordingly.

In addition, a novel method is presented that is basically along the lines of Hartree et al. (1946) but employs a novel differential equation, in unapproximated form. This method is found to be robust, easy to implement without any pitfalls, and accurate under the most general conditions. To this end, the new method is compared to the ⅘ earth radius model as well as to our revised version of the method of Caumont et al. (2006), representative for the class of methods employing Snell’s law.

During the writing of this manuscript, the work of van der Werf (2003) came to the attention of the authors, which derived and solved yet another system of differential equations.

First, we recall the basics of radar beam propagation conditions in section 2, with special attention given to so-called ducts in section 3. Sections 4 and 5 describe the above-mentioned two methods that we chose for comparison, and section 6 introduces the new method. In section 7, sensitivity experiments are carried out and the results are analyzed.

## 2. Radar beam propagation

A radar beam that is propagated through the atmosphere encounters variations of refractive index along its trajectory, which cause the beam to become curved. The total angular refraction of the beam between two points is commonly called bending. Here, what we call “beam” is understood as a collection of neighboring rays that individually may undergo different refraction, leading to a distortion of the beam.

It is helpful to briefly recall the physical basis for computation of atmospheric refraction at first. More thorough treatments on the subject may be found in, for example, Kerr (1951) and Bean and Dutton (1966), as well as Doviak and Zrnić (1993).

For describing the ray path, the classical geometric optics is a commonly used approximation. This approximation is applicable, if, within one wavelength of radiation,

the refractive index

*n*changes only very little andthe mutual distances between neighboring rays change also only very little.

Under these conditions, a single ray path is determined by Fermat’s principle, which states that the travel time *t* between two points *A* and *B* be minimal. Travel time depends on the propagation speed, which is given by *c*′ = *c*/*n*, where *c* is the speed of light in vacuum and *n* is the refractive index of the medium (to be precise: its real part).

Fermat’s principle reads

where *dr* is an infinitesimal arc element. This is a classical problem of functional analysis. For the atmospheric ray propagation, it is assumed that the earth is spherical with radius *R*_{E} [defining the above mean sea level (MSL)]. The refractive index *n*, in general, varies in all spatial directions, but in the atmosphere, vertical variations are usually much larger than horizontal variations. Therefore, it is assumed hereafter that *n* only depends on height *h* over the earth’s surface.

Under these conditions, the ray path is defined by the Euler–Lagrange equation of the system, which reads that after transforming *dr* to the arc distance element *ds* at MSL height *h* = 0 (Hartree et al. 1946),

It is worth noting that Eq. (2) does not allow for exact vertical propagation because in that case *dh*/*ds* → ∞. Moreover, one can show that this second-order nonlinear ordinary differential equation (ODE) is almost equivalent to the following integral conserved quantity:

where is the local elevation and is given by

By “almost” we mean that integration is mathematically not an equivalent transformation and additional (nonphysical) solutions can be created by integration. Here, this is manifested by a sign ambiguity of *ϵ* in Eq. (3), because cos*ϵ* = cos(−*ϵ*). We will address this problem later in section 5.

Equation (3) is the well-known Snell’s law for a continuous spherically stratified medium, which states that the constant on the rhs is conserved along a ray path. This conserved quantity has a similar significance as, for example, the mechanical energy for the equation of motion.

A useful simplified approximation of Eq. (2) may be obtained by assuming: *R*_{E} + *h* ≈ *R*_{E}, *dh*/*ds* ≪ 1 (rays at low elevations), and *n* ≈ 1, so that Eq. (2) becomes

Based on Eq. (5), Doviak and Zrnić (1993) showed that the curvature of the ray *C*_{0} is

Since *dh*/*ds* ≈ *ϵ* for small *ϵ*, the term *d*^{2}*h*/*ds*^{2} describes the change of the local elevation with *s*. Hence, it is clear from Eq. (5) that if *dn*/*dh* = −1/*R*_{E}, then *dh*/*ds* is constant and equal to zero if *ϵ* = 0; that is, the ray spreads parallel to the earth’s surface and the curvature of the ray path is 1/*R*_{E}, but its curvature relative to the earth is zero. With Eq. (6) we can conclude that the ray’s curvature *C*_{E} relative to the earth is

## 3. Refractive index

The variation of *n* is closely related to the vertical variation of temperature *T*, water vapor partial pressure *e*, and total pressure *p*. As *n* is slightly larger than unity (e.g., 1.0003), it is much more convenient to define the so-called refractivity *N*, given by

For instance, if *n* is 1.0003, then the corresponding value of *N* is 300. Bean and Dutton (1966) showed that *N* can be empirically given by

where *c*_{1} = 77.6 K hPa^{−1}, *c*_{2} = −6.0 K hPa^{−1}, *c*_{3} = 3.75 × 10^{5} K^{2} hPa^{−1}.

Normally, *N* decreases with height in the atmosphere, which leads to a downward bending of radar beams. Under some circumstances, *N* may increase with the height—that is, *dN*/*dh* > 0—and the beam bends away from the earth’s surface. As mentioned above, a horizontal ray (*dh*/*ds* = 0) remains horizontal (i.e., has the same curvature as the earth), if *dn*/*dh* = −1/*R*_{E} respectively *dN*/*dh* = −10^{6}/*R*_{E}, and a nonhorizontal ray would preserve its local elevation. For example, if *dh*/*ds* > 0, then a quasi-helical motion around the earth would result. With *R*_{E} ≈ 6371 km, this is *dN*/*dh* = −157 km^{−1}. If the derivative of *N* is smaller than that, then the curvature of the ray becomes larger than that of the earth, and the ray will, after reaching a maximum height, be bent down and trapped between this height and the earth’s surface. This process is called trapping, and the layer of the atmosphere within which the beam is bent back downward is called a trapping layer. If there is a region below the trapping layer with a larger derivative of *N*, then the mode of beam propagation is similar to that of a waveguide, and this configuration is called a duct.

For the lowest 1 km of the International Civil Aviation Organization (ICAO) standard atmosphere, Doviak and Zrnić (1993) give a value of *dN*/*dh* = −40 km^{−1}, and we will refer to this as normal conditions in the following, as they represent some “average” climatological conditions near the ground.

If *dN*/*dh* lies between 0 and −40 km^{−1}, then the beam will be bent toward the earth’s surface with a curvature less than that of the normal conditions, and we refer to it as *subrefraction* (see Fig. 1). Superrefraction occurs when *dN*/*dh* ranges from −40 to −157 km^{−1}. In this situation the beam is bent down to the surface at a rate less than the earth’s curvature but more than normal.

When considering atmospheric ducts, instead of *N* the so-called modified refractivity *M* is usually preferred, defined as

Then *dM*/*dh* = *dN*/*dh* + 157 km^{−1} and for constant *M* (*dM*/*dh* = 0) the curvature of the propagation of a nearly horizontal beam is that of the earth’s surface, and *dM*/*dh* < 0 for trapping conditions. Figure 1 shows the various categories of refraction in terms of *dN*/*dh* and *dM*/*dh*.

According to the profile of *M*, three basic forms of a duct with corresponding duct depths are shown in Fig. 2. The case in Fig. 2a illustrates the structure associated with a simple surface duct. Here, the duct extends from the local minimum to the surface, and the trapping layer, where *dM*/*dh* < 0, stretches throughout the duct. Figure 2b is referred to as the surface S-shaped duct, which reaches down to the surface, while the trapping layer does not, since *dM*/*dh* > 0 near the surface. In these two cases, the duct depth is the height difference between the ground and the top of the duct, where the minimum in modified refractivity profile is achieved. In Fig. 2c, the common conditions for an elevated duct are given, where the value of *M* at the surface is less than that at the top of the duct, and so the duct cannot reach down to the surface. Its depth extends from the local minimum to the height at which the *M* value equals that at the top of the duct.

As mentioned above, a duct is the result of strong vertical changes in the refractive index of the atmosphere between air masses of different temperatures and humidities, especially at low levels of the atmosphere. A duct can occur on a very large scale when a large mass of cold air is overrun by warm air, leading to a strong temperature inversion. A duct can also occur when a strong cap of warm and dry air exists in the lower troposphere above very moist air. On one hand, a duct causes the electromagnetic energy to be able to propagate over farther distances, allowing long-range radio communication; on the other hand, in weather radar applications, ducts usually lead to coverage fades, increased ground clutter, increased anomalous propagation, and range–height errors. One part of the radar simulation process is the computation of beam propagation within the atmosphere simulated by an NWP model in an appropriate way. It is known that low elevations are often vulnerable to anomalous propagation and orographic beam blockage, which can seriously affect the radar’s ability to detect and quantify precipitation at ground level. Important issues here are to minimize the influences of these effects in the observation–simulation comparison. In the following sections, we briefly describe and analyze a simple well-known approximate technique and two more sophisticated (new) methods.

## 4. 4/3 earth radius model

For convenience of computation, in this model the geometry is transformed in such a way that the earth has a fictive radius/curvature and the ray propagates straightforward. However, the relative curvature of the ray toward the earth is retained; that is

which yields

Here, *R*_{eff} represents the effective earth’s curvature relative to a straight ray and *K*_{eff} is the effective earth radius factor depending on *dn*/*dh*. From Eq. (12) we can see that if *dn*/*dh* is constant, then the earth has an effective radius of constant *R*_{eff}.

In the standard atmosphere, where the refractive index decreases linearly with the height in the mean by *dn*/*dh* = −40 × 10^{−6} km^{−1} in the lowest 1 km or so (i.e., −40 *N*-units km^{−1} or 117 *M*-units km^{−1}), then we obtain for *K*_{eff},

This is a common model to approximate ray paths, which assumes that the effective earth radius is one-third larger than the real one, so *R*_{eff} = 4/3*R*_{E}. This model allows for a straightforward analytical estimation of each pulse volume height *h* and surface distance *s* relative to the radar site at a height of *h* = 0 (MSL) (cf. Fig. 3), given as

This model is referred to as the 4/3 earth radius model and is abbreviated in the following as 43ERM. As shown in Doviak and Zrnić (1993), for weather radar applications 43ERM can be used for all elevations if *h* is confined to the first 10–20 km and if *n* has a slope of −1/(4*R*_{E}) in the first kilometer of the atmosphere. But the slope of *n* is usually not constant. If *n* decreases much more rapidly than in the standard atmosphere, then the beam will likely be bent downward, and then the height of pulse volumes tends to be overestimated by 43ERM. These errors can be quite significant for elevations smaller than ~1°. More sophisticated methods, taking into account the actual atmospheric conditions instead of climatology, can be found in Caumont et al. (2006), Gao et al. (2006), Chen et al. (2009), and van der Werf (2003). Caumont et al. (2006) represents a class of methods based on Snell’s law, Eq. (3), that is reviewed and revised in the next section.

## 5. Method based on total reflection

Under realistic atmospheric conditions—for example, an *n* profile based on a radiosonde measurement—some authors computed the ray propagation iteratively by discretizing Snell’s law, Eq. (3), in the along-beam direction *r* in steps of some fixed increment Δ*r*.

In this sense, a method used in Caumont et al. (2006) has been revised, adding suitable criteria for catching inflection points under ducting conditions (“total reflection”) and for negative antenna elevations under nonducting conditions (the effect of earth’s curvature).

Let *l* = 1, …,*L* be the numeration index of steps. Then the height *h*_{l} and the MSL-reduced surface distance *s*_{l} at some location are iteratively calculated from the values at *l* − 1 (height above MSL *h*_{l−1}, local elevation *ϵ*_{l−1}) under the assumption of straight rays within a prespecified range step Δ*r*, written as

Figure 4 shows a sketch of these quantities.

Given the new *h*_{l}, the new *ϵ*_{l} can be directly derived from the discretized Eq. (3) using

where *ϵ*_{l} is the local elevation of the ray at range *l*Δ*r*.

Before this discretized iterative method is further developed below, it is instructive to highlight an important property of the “true” ray path—the analytic solution, if you will—concerning inflection points under ducting conditions. From Eq. (3) and assuming that *n*(*h*) is continuous, we can derive an analytic implicit equation for the height *h*_{*} at which *ϵ*_{*} = 0 (inflection), whose lhs resembles *F* in Eq. (18),

where the index 0 denotes the starting values at the antenna. Note that the nominator on the lhs has exactly the same value as the one in *F*. If this equation has at least one real solution for *h*_{*}, then the ray path exhibits at least one inflection point at *h*_{*}. If *n*(*h*) is a linear function *n*(*h*) = *n*_{0} + *γ*(*h* − *h*_{0}), where *γ* is constant, then there can be one inflection point at the most. Inserting *n*() = *n*_{0} + *γ*( − *h*_{0}) into Eq. (19) leads to a quadratic equation for *h*_{*} that can be solved explicitly,

with

The subscript “1/2” indicates two possible solutions. If the roots are real numbers, then the inflection point *h*_{*} is one of them.

From this it follows that the *F* in equation Eq. (18) becomes 1 if step *l* hits exactly an inflection point, and in the limit Δ*r* → 0, and that if the ray path is continuous, *F* cannot become larger than 1. However, for finite Δ*r* the *F* can become larger than 1, which is a discretization artifact.

Now let us turn back to Eq. (18). Two problems arise here:

arccos is not defined in the case of

*F*> 1 andthe sign of

*ϵ*_{l}is ambiguous (±) because arccos is not a unique mapping for the codomain [−*π*/2,*π*/2].

Concerning problem 1, this could happen in our discretiztation if *n*_{l} is “sufficiently” smaller than *n*_{l−1} at some location *l*. In analogy to total reflection at a discrete *n*-jump, we here assume that the ray exhibits an inflection point and is reflected back internally, so

Term Δ*ϵ* is a correction term that accounts for the effect of the earth’s curvature on the local elevation along a Δ*r* segment. A graphical derivation of Δ*ϵ* can be found in Fig. 5.

If *F* < 1, then it is reasonable to assume that the sign of *ϵ* does not change from one step to the next, so we have instead of Eq. (18) the following:

Note that we include the case *F* = 1 into the forced sign change, although the arccos would be defined. This is to avoid an artificial numerical trapping of the ray at *h*_{*} under ducting conditions, which can be identified by explicitly expanding the iteration algorithm equations (16) and (22) but replacing the condition *F* < 1 in Eq. (22) with *F* ≤ 1.

Next, we give an insight into when *F* ≥ 1 can occur in the discretized equations. From the physical understanding, total reflection occurs when *ϵ*_{l−1} is very small, meaning *ϵ*_{l−1} ≈ 0 and cos*ϵ*_{l−1} ≈ 1, thus

which means *F* ≥ 1 when Δ*n* + Δ*h*/*R*_{E} ≤ 0 or

holds. For the case of positive elevations, which means the ray propagates upward, we have Δ*h* > 0 and can rewrite Eq. (24) as

denoting ducting conditions (see Fig. 1). Equation (25) says that for the ascending ray, *F* ≥ 1 is possible if *ϵ*_{l−1} ≳ 0 and *dn*/*dh* ≤ −157 × 10^{−6} km^{−1}. It is analogous for the case of negative elevations: the ray goes downward, Δ*h* < 0, and *F* ≥ 1 is possible when *ϵ*_{l−1} ≲ 0 and *dn*/*dh* ≥ −157 × 10^{−6} km^{−1}, that is, nonducting conditions.

Note, however, that some approximations are involved to arrive at Eq. (24). In practice, we found by extensive testing the appearance *F* ≥ 1 under ducting conditions is only robust for positive elevations. Conversely, it almost never happens for negative elevations under nonducting conditions, and the sign change from negative to positive because of the earth’s curvature fails to occur. As a consequence, *ϵ*_{l} gets trapped at very small negative angles and the beam remains essentially parallel to the earth’s surface (spurious ducting).

This behavior is investigated more closely in the appendix. It is identified as an artifact of the forward discretization and an ad hoc criterion to catch it is motivated. The criterion uses the increment between *ϵ*_{l−1} and *ϵ*_{l−3} to linearly extrapolate and predict *ϵ*_{l}. A sign change is assumed if *ϵ*_{l−1} < 0 and *ϵ*_{l−1} + (*ϵ*_{l−1} − *ϵ*_{l−3}) > 0. This criterion has been tested for many different *n* gradients and negative antenna elevations, and is also beneficial in experiment 4 in section 7.

Therefore, Eq. (22) is modified as follows:

Because of the newly considered total reflection assumption, this modified method is called TORE (acronym for total reflection) and can be summarized as follows (see Fig. 6):

Step 1: Calculate height

*h*_{l}and MSL-reduced surface distance*s*_{l}using Eqs. (16) and (17), starting from*h*_{l−1}and*s*_{l−1}with local elevation*ϵ*_{l−1}.Step 2: Estimate

*n*_{l}=*n*(*h*_{l}) using radiosonde or NWP data.- Step 3: Calculate
*ϵ*_{l}using Eq. (26). Note that*n*_{l−1}is known from step 2 of the previous iteration.Fig. 6.

Steps 1–3 are repeated from *l* = 1 to *l* = *L*. In the first iteration, the values at *l* − 1 are antenna elevation *ϵ*_{0}, height *h*_{0} at the radar antenna, refractive index *n*_{0} = *n*(*h*_{0}), and *s*_{0} = 0.

Note that despite extensive testing, we cannot exclude that the above-mentioned ad hoc criterion might fail in rare instances because it is not rigorously mathematically well founded. Note also that the above-mentioned sign ambiguity is a general problem of Snell’s law (as stated earlier), and that all methods based on it have to deal with the problem in one way or another.

## 6. Method using a Second-Order Ordinary Differential Equation

Although TORE considers explicitly the actual refractive index, an ad hoc criterion is required to determine the sign change of local elevations, and we cannot guarantee that this criterion is robust under all possible conditions. In this paper, we have found a novel method, called SODE, that offers a straightforward analytical–numerical solution for ray propagation and considers the sign change automatically.

As can be seen from Eq. (2), the ray propagation can be formulated as an initial value problem (IVP) of an ODE. Principally this would be possible by employing Eq. (2) directly, but it has the drawback of being formulated relative to the independent coordinate *s*. For many practical applications—for example, as part of a radar forward operator—a formulation relative to the along-beam range *r* would be preferred because it is then possible to discretize the solution by using a constant Δ*r* that is directly related to the radar range. To derive such an alternative ODE, we start from Snell’s law, Eq. (3), but now assume *h* as a function of *r*, that is,

From differential geometry (cf. Fig. 7) one obtains for infinitely small *dh* and *dr*

therefore,

As indicated, *h* is assumed to be a function of the range *r*, so the refractive index *n* depends implicitly on *r*. One differentiates Eq. (30) with respect to *r* and obtains

The singularities *dh*/*dr* = ±1 can be removed by multiplying , which yields

Since *dh*/*dr* ≡ 0 as a general solution is unphysical, it holds

One divides both sides of Eq. (33) by *n*(*R*_{E} + *h*) and obtains the second-order nonlinear ODE

and, by substituting *dh*/*dr* = *u*, the equivalent set of two coupled first-order equations are arrived at

Equation (34) and the equation system (35) are physically equivalent to Eq. (2), but are formulated with independent coordinate *r* instead of *s*. It can be easily checked that Eq. (34) also includes the constant solution *dh*/*dr* ≡ 0 describing a ray along the earth’s surface for *dn*/*dh* = −*n*/(*R*_{E} + *h*) and the cases *dh*/*dr* = ±1 for the exact vertical propagation [an advantage against Eq. (2)].

Equation (35) can be treated as an IVP with initial values

and the ray tracing problem is then uniquely solved by this IVP.

In analogy to TORE, Eq. (35) is discretized and solved in steps of Δ*r*. The iteration step from location *l* − 1 to *l* is performed as follows:

Step 1: Estimate l/

*n*_{l−1}and*dn*/*dh*|_{l−1}at height*h*_{l−1}using radiosonde soundings or NWP data.Step 2: Solve Eq. (35) with initial values

*u*_{l−1}and*h*_{l−1}to obtain*u*_{l}and*h*_{l}*.*- Step 3: As in Eq. (17), calculate MSL-reduced surface distance
*s*_{l}from with Steps 1–3 are repeated from*l*= 1 to*l*=*L*. Note that the IVP posed in step 2 is currently solved using the fourth-order explicit Runge–Kutta method (hereafter RK4), but any other numerical standard method for solving ODEs would be suitable as well. For the first iteration*l*= 1, the initial values in step 2 are given using Eqs. (36) and (37), and the refractive index and its slope have to be estimated at*h*_{0}, resulting in 1/*n*_{0}and*dn*/*dh*|_{0}, respectively.

Other authors have also applied differential equation solvers for the ray propagation problem (Hartree et al. 1946; van der Werf 2003), but the formulation of the above-mentioned ODE in terms of *r* is believed to be new and especially suitable for radar forward operators.

## 7. Sensitivity experiments

So far we have presented three methods to calculate radar beam propagation in a stratified atmosphere. In what follows, we compare all these methods by evaluating them for specific atmospheric conditions. This is done with a series of sensitivity experiments in a framework, where certain horizontally homogeneous vertical profiles of *T*, *p*, and *e* are prescribed (M. Neuper 2010, unpublished manuscript). The first three experiments are based on the idealized ducting profiles introduced in Fig. 2. A fourth experiment is based on standard atmosphere data, and a fifth experiment applies measured radiosonde data from a ducting case.

For all experiments we choose a maximal surface cover range of 300 km and a range resolution Δ*r* of 500 m. In the first three experiments, we investigate simulations of beam propagation for two initial elevations, *ϵ*_{0} = 0.1° and *ϵ*_{0} = 1.1°, and in the fourth one for *ϵ*_{0} = −0.3°. To stimulate different kinds of ducts, the radar antenna is set accordingly to different heights in the experiments.

### a. Experiment 1: Idealized surface duct

In this experimental setup, we simulate a surface duct. Accordingly, we have chosen the profiles of *M* and *N* with respect to *h* as given in Figs. 8a and 8b, respectively. The radar antenna is chosen at a height of 200 m. Figure 8a shows a large negative slope of −100 *M*-units km^{−1} of *M* for the lowest 350 m and thereafter a slope of 117 *M*-units km^{−1}.

The simulation results are shown in Figs. 8c and 8d. Figure 8c represents the variations of beam heights *h* with distance *s* computed using the three methods, while Fig. 8d shows the absolute height differences of TORE and SODE compared to 43ERM. One can see that for *ϵ*_{0} = 1.1°, the beam heights calculated using all three methods are generally close to each other (as shown in Fig. 8d) and differ less than 600 m at maximum range; for *ϵ*_{0} = 0.1°, the resulting rays according to TORE and SODE are refracted downward to the earth’s surface because of the strong negative *dM*/*dh*, while 43ERM produces a curve that is straightens upward No surface reflection was taken into account in these calculations, and therefore the rays of TORE and SODE end when they reach the ground. The discrepancies compared to 43ERM grow already to about 1600 m at a distance of 140 km. This observation is therefore in accordance with the statement from Doviak and Zrnić (1993) that for the higher elevations the radar beams are less sensitive to the refractivity gradient, while for low elevations (<1°) under ducting conditions, 43ERM is prone to (strongly) overestimating the beam heights.

### b. Experiment 2: Idealized surface S-shaped duct

Now we apply all these methods to an idealized surface S-shaped duct, characterized by the profiles of *M* and *N*, shown in Figs. 9a and 9b. The profile of *M* begins with a slope of 117 *M*-units km^{−1} for the lowest 100 m and then alters to −100 *M*-units km^{−1} until 400 m, after which it goes back to 117 *M*-units km^{−1}. The antenna’s height is chosen to 40 m, slightly below the trapping layer. As can be seen in Figs. 9c and 9d, for *ϵ*_{0} = 1.1° the differences of the three methods are bounded to 1 km at the maximum distance. However, for *ϵ*_{0} = 0.1°, the rays calculated by both SODE and TORE propagate in a wavelike mode, trapped within the ducting layer (see Fig. 9c), and the height differences compared to 43ERM reach already 1600 m at a distance of 150 km (see Fig. 9d). The reason for the slight discrepancies between SODE and TORE for *ϵ*_{0} = 0.1° in Fig. 9c will be addressed in experiment 5.

### c. Experiment 3: Idealized elevated duct

A further situation modeled is an idealized elevated duct. The corresponding profiles of *M* and *N* are illustrated in Figs. 10a and 10b, respectively. The profile of *M* starts with a slope of 117 *M*-units km^{−1} for the first 250 m and then changes to −100 *M*-units km^{−1} until 400 m, and at last returns to 117 *M*-units km^{−1}. The antenna height is 300 m within the trapping layer. The general features of the results, illustrated in Figs. 10c and 10d, are mainly the same as those in experiment 2, except that the duct here is lifted in the air and does not touch the ground.

### d. Experiment 4: Standard conditions

This experiment is based on data for standard conditions. As shown in Figs. 11a and 11b, now a constant slope of *dM*/*dh* = 117 *M*-units km^{−1} throughout the atmosphere is considered. The antenna height is set to 200 m. But now the antenna elevation is set to a negative value, *ϵ*_{0} = −0.3°. To demonstrate the effects of the ad hoc criterion in TORE for sign change [cf. Eq. (26)], two experiments are performed here: one using Eq. (26), denoted with E4(1), and the other one using Eq. (22) without this criterion, denoted with E4(2). In E4(1), the results of all three methods are basically identical (see Figs. 11c and 11e). All prompt the beams to descend at the beginning due to the negative initial elevation and to slope upward after a distance of about 50 km due to the earth’s curvature. This shows that, for normal conditions, 43ERM is a satisfactory approximation in comparison with the solution of SODE, which is considered as an accurate reference solution.

But in E4(2), TORE is not able to overcome the slight negative elevations (near to 0° as shown in Fig. 11e) and thus flattens out afterward. The reason is explained in the appendix and is basically an artifact of the forward discretization of the ray path under nonducting conditions and for negative elevation angles. Here, the tangentially discrete ray segment defined by Eq. (16) shoots to the “wrong” side of the convex ray path and does not reach a point where *F* ≥ 1. Thus, the conditions for a sign change are not met and the negative local elevations fail to become positive because Eq. (22) preserves the sign. This shows the necessity of the ad hoc criterion in Eq. (26) as a further check for the sign change of ∈ from negative to positive just because of the earth’s curvature.

Note that experiments 2 and 3 (quasi-standard conditions in the lowest levels) have been also tested without the ad hoc criterion, and that the TORE rays have also gotten trapped horizontally at and after their first height minimum.

### e. Experiment 5 using measured radiosonde data

As a last experiment, a case based on real atmospheric conditions is investigated, one that exhibits a strong temperature inversion and moisture profile near the ground observed at Stuttgart–Schnarrenberg (WMO station 10739) in Germany at 0000 UTC 4 September 2004 (see Fig. 12). The corresponding profiles of *M* and *N* given in Figs. 13a and 13b, respectively, are derived from the radiosonde data available from the University of Wyoming,^{1} where a vertical interpolation of original *T*, dewpoint *T*_{d}, and *p* data to additional levels (linear oversampling every 10 m vertical) is performed and from these oversampled data, *n* is computed. The refractive index *n*_{l} at some arbitrary level *l* is derived by linear interpolation from upper and lower neighboring oversampling points at locations *l*_{>} and *l*_{<}, respectively. The refractivity slope is approximated by a simple differential quotient,

The linear oversampling of the original radiosonde data minimizes interpolation artifacts for *n* and its vertical slope. It is justified because in radiosonde data, *T* and *T*_{d} often exhibit a near-linear dependence on height in between the data points, which are stored in radiosonde datasets, and *p* varies smoothly with height.

If, however, some noise in the derived *n* profiles should lead to noisy gradients, then some smoothing could be obtained by applying more sophisticated methods for interpolation and slope calculation. The authors found the so-called Savitzky–Golay filter (Press et al. 1992) very useful, that is, fitting of a low-order polynomial to a wider stencil of neighboring oversampling points (e.g., fifth-order polynomial, 12 surrounding points) and computing *n* and its slope from this polynomial instead of interpolating from the original data.

The setup of maximum surface distance, range resolution, and initial elevations remains the same as previously and the radar antenna height is set to 40 m.

The corresponding *M* profile in Fig. 13a shows a duct between 17 and 144 m above the surface (i.e., duct depth = 127 m) and a trapping layer extending from 71 to 144 m. Note the thin layer with a positive *N* gradient below the duct. Therefore, it can be expected to observe an elevated duct between 17 and 144 m in ray paths with low elevations. Figure 13c illustrates the comparison of beam heights computed by the three methods. As expected, all three methods provide nearly the same results for *ϵ*_{0} = 1.1°; for *ϵ*_{0} = 0.1°, 43ERM, not “knowing” about the actual ducting conditions, generates a lifting curve, while TORE and SODE are consistent with the expected ducting conditions, and both are able to deliver reasonable waveguide-like results but with slightly different shapes (Figs. 13c and 13d). Figure 13e shows the corresponding local elevations.

The discrepancies between SODE and TORE seen in Fig. 13c arise from the different trajectories of *ϵ*_{l} obtained by both methods (see Fig. 13e). Overall, TORE simply changes the sign of *ϵ*_{l} at inflection points, while the change of *ϵ*_{l} by SODE is more continuous and thus slower because of the differentiability of solution. Some of the changes in *ϵ*_{l} for SODE seem to be abrupt because the corresponding solid curve is plotted only at each twentieth point of the original solution. Nevertheless, SODE still results in a periodic curve for *ϵ*_{l} as TORE does but with a longer wavelength, and both theoretically converge to the same limit as the resolution Δ*r* approaches zero. Therefore, we assume that the aforementioned discrepancies could be related to some stability issue of the numerical method (RK4) used for solving SODE and/or to the resolution Δ*r*. As a matter of fact, these discrepancies almost disappear if we replace Δ*r* = 500 m with Δ*r* = 200 m, as illustrated in Figs. 14c and 14e. For a further insight into the convergence behavior in terms of Δ*r*, more experiments in the same framework have been done with Δ*r* ranging from 2 to 4 km (not shown here). It has been found that TORE with the ad hoc criterion still works for the very small step size of Δ*r* = 2 m. Toward very large step sizes, SODE fails to produce a ray path in a waveguide-like fashion beginning at Δ*r* = 2 km and TORE beginning at Δ*r* = 4 km, which may also indicate a better stability of TORE in this case. However, we would like to point out that the convergence behavior with varying Δ*r* is very dependent on the *n* profile; therefore, the range of Δ*r* for convergence may vary from case to case. Regarding this, an analytic formula that links the convergence range of Δ*r* to properties of the *n* profile as well as *ϵ*_{0} would be desirable but is beyond the scope of this paper.

It is worth noting that TORE without the ad hoc criterion also works correctly in this case and does not produce a spuriously trapped ray after the first height minimum (not shown here). This, however, is only the result of a “lucky” configuration, namely, the height layer with increasing *N* until 30 m, which is sufficient for the forward iteration to reach a point where *F* ≥ 1, enabling the local elevation to change its sign.

## 8. Summary and discussion

In this paper, we assessed the performance of three radar beam tracing methods by several sensitivity experiments: 43ERM (well-known and based on atmospheric standard conditions); TORE (partly known from literature and based on actual vertical profiles of refractive index *n*); and SODE (new method, also based on actual profiles of *n* and introduced in section 6). Both the TORE and SODE methods employ actual *n* data and are rigorously based on geometrical optics and its fundamental Fermat’s principle. Whereas SODE involves the solution of a novel ordinary differential equation in the form of an initial value problem, TORE is based on the conservation of an integral quantity of this differential equation along the ray path, known as Snell’s law for continuously stratified media.

Because 43ERM does not take into account the true environmental conditions at a specific time, it tends to overestimate the beam heights in the case of superrefraction or ducting, especially for low elevations (*ϵ*_{0} ≲ 1.0°). However, for conditions that are near normal in the lowest 1000 m or so of the atmosphere or for higher antenna elevations, it generally works well, as also noted in earlier studies (e.g., Doviak and Zrnić 1993) and is a commonly used method among radar specialists today. Such conditions are prevailing for the vast majority of cases because most radar data are taken from elevations > 1° and because superrefraction or ducting conditions occur only occasionally.

However, when it comes to superrefraction or ducting conditions connected with low elevations, TORE is able to reproduce the ducting effect. In order for this to work, the original method of Caumont et al. (2006) had to be extended—in light of the ambiguity of the arccos function regarding the sign of the local elevation angle—by a proper method and criteria to detect and treat inflection points correctly, where the sign changes. This happens either for upward-propagating rays in ducting conditions or downward-propagating rays in nonducting conditions. For the latter case, our new ad hoc approach to account for elevation sign changes from negative to positive plays a key role, without which the beam would get artificially trapped at a local elevation *ϵ* ≲ 0° and would propagate purely horizontally farther out. Nevertheless, this ad hoc approach is based on very simple linear extrapolation, and its stability and robustness, although satisfactory in all our cases, have not been rigorously mathematically proven and cannot be guaranteed.

Instead, inflection points occur naturally with the new SODE method, which means that no special measures are necessary to correctly detect and treat them. And it also provides reasonable and robust results in all presented tests. Compared to the traditional beam propagation, Eq. (2), SODE is also preferable, for example, in the implementation of radar forward operators because of its dependency on *r* instead of *s*.

It can be said that SODE is more robust than TORE if accounting for the actual atmospheric conditions because of its straightforward numerical solution and automatic consideration of inflection points. Therefore, we regard SODE as a reference method. However, prerequisite for a successful detailed beam propagation computation is a very good knowledge of the 3D atmospheric state—that is, *n* and its vertical slope—which may vary also horizontally (this last point is not taken into account in the present paper). With today’s aerological network (sparse number of stations and sparse observation times), this is certainly not the case in general, and the results can only be as good as the input data. However, if one day better—that is, spatially and temporally more dense—observations should be available, then SODE can play out its advantages. But already today it could be useful—owing to its robustness—in automated software tools for propagation computations, as well as for radar forward operators, which compute synthetic radar data in the virtual world of the simulation results of atmospheric simulation models. Concerning the last point, most (operational) NWP models today lack the necessary vertical resolution for accurate beam propagation computation, but for the future it is foreseeable that much effort will go into higher model resolution (at least in research), so that then SODE could be the method of choice.

## Acknowledgments

The authors thank the German Weather Service (DWD) for funding this work under its Extramural Research Program. Thanks also to Prof. Dr. Klaus D. Beheng for his critical comments on the first draft, and to three anonymous reviewers for their valuable corrections. Especially one anonymous reviewer provided very constructive input and pointed us to a problem in our first implementation of the TORE method. This motivated us to find the ad hoc solution mentioned in section 5.

### APPENDIX

#### Spurious Ray Trapping in TORE for Negative Antenna Elevations

The problem of spurious ray trapping in the TORE method for negative elevations under nonducting conditions, which has been mentioned in sections 5 and 7 (experiment 4), is further investigated here. We specifically refer to the atmospheric conditions of experiment 4, subexperiment E4(2), of section 7.

Under these conditions, the height *h*_{*} in Eq. (19), which denotes the true height of the inflection point, exists and can be explicitly computed from Eq. (20) because *n*(*h*) is linear. During the TORE iteration, *h*_{l} generated by Eq. (16) approaches *h*_{*} from the positive side as *ϵ*_{l−1} approaches 0 from the negative side in a polygonal chain with *h*_{l−1} > *h*_{l} > *h*_{*}, up to the point where *ϵ*_{l−1} becomes sufficiently small to principally enable *F* to become equal or larger to 1, but it does not actually happen. This is depicted in Fig. A1. In contrast to experiment 4, we chose *ϵ*_{0} = −0.1° instead of −0.3°. The top panels show the discrete *ϵ*_{l} as function of surface distance, as crosses and diamonds, zoomed to the location where the sign change is to be expected, for Δ*r* = 400 m (left, solid curve) and Δ*r* = 200 m (right, dashed curve), respectively. The bottom panels depict the actual polygonal chains of these two discretized rays in true Cartesian coordinates, which are referenced in a vertical cut along the ray path and through the (spherical) earth’s center:

with *h*′ and *ϕ*′ chosen in such a way that the coordinate center is around the expected inflection point. In this way, the segments of the discrete polygonal ray paths are geometrically undistorted with respect to true straight rays, as one would see it from space. This will be helpful for the following discussion. The left panel shows only the polygon for Δ*r* = 400 m, whereas the right panel contains both. The gray line denotes *h*_{*} (*F* = 1) and is parallel to the earth’s surface. For points above this line we have *F* < 1 and below we have *F* > 1. Note that the *y* axis is extremely stretched and spans only ±2 cm, whereas *x* spans ±600 m. Elevation angles at the crosses and diamonds have to be measured with respect to parallels to the gray line and are all negative, corresponding to the graphs in the top panels.

To trigger a sign change, the *h*_{l} (crosses) would have to hit the gray line at some point or even fall below. In the bottom-left panel, we can see why this is not happening. Consider the step from *l* − 1 to *l*. The true (analytical) curve would be something like the convex dashed line, which touches the gray line at some location and causes the sign change. But in the discretization, this convex line is approximated by its tangent at *l* − 1 and *h*_{l} overshoots *h*_{*}. Therefore, no sign change occurs; the next step proceeds with a small negative elevation but again overshoots *h*_{*} and so on. The ray polygon is trapped parallel to the gray line and is parallel to the earth’s surface.

At half Δ*r* (dashed curve, right panels), the problem basically remains the same, although the absolute value of the trapped elevation angle is about half. Further reducing Δ*r* does also not solve the problem, only that the trapped elevation proportionally approaches 0. The effect occurs also for any other nonducting *n* gradient (for ducting conditions, there is no sign change in case of negative antenna elevations).

All in all, the problem can be identified as a forward discretization issue. If the region where *F* ≥ 1 is toward the inside of a convex ray path, then the forward tangential approximation shoots to the other side and does not hit a point in space with *F* ≥ 1. Conversely, for upward-propagating rays and ducting conditions (region with *F* ≥ 1 toward the outside), the tangential shoots to the correct side and can hit a point with *F* ≥ 1 and the sign can change.

To circumvent this problem, an extra ad hoc criterion is proposed. Because the absolute value of the trapped elevation scales with Δ*r*, it seems reasonable to predict a sign change based on a simple linear extrapolation of the previous elevation changes. But we can see in the top-left panel that the increment between *ϵ*_{l−1} and *ϵ*_{l−2} is not enough to get a positive *ϵ*_{l} = *ϵ*_{l−1} + (*ϵ*_{l−1} − *ϵ*_{l−2}). However, replacing *ϵ*_{l−2} with *ϵ*_{l−3}, as marked by dashed lines in the figure, does the trick.

Therefore, a sign change is assumed if *ϵ*_{l−1} < 0 and *ϵ*_{l−1} + (*ϵ*_{l−1} − *ϵ*_{l−3}) > 0. This criterion has been successfully tested for various negative elevations and *n* gradients.

A more general solution could be to formulate the TORE iteration in an implicit way instead of explicit, but this would likely complicate the computations.

## REFERENCES

*Radio Meteorology*. National Bureau of Standards Monograph, No. 92, U.S. Government Printing Office, 435 pp.

*Doppler Radar and Weather Observations*. 2nd ed. Academic Press, Inc., 562 pp.

*Propagation of Short Radio Waves,*D. E. Kerr, Ed., McGraw-Hill, 41–58.

*Meteorological Factors in Radio-Wave Propagation,*The Physical Society, 127–168.

*Numerical Recipes in Fortran 77: The Art of Scientific Computing,*Vol. 1,

*Fortran Numerical Recipes,*Cambridge University Press, 933 pp.

## Footnotes

Current affiliation: Deutscher Wetterdienst, Offenbach, Germany.

^{1}

Although a radiosonde measures on a spatial scale much larger than the radar wavelength, the data are readily available and are a good source of atmospheric information on temperature and humidity structure. It is assumed that the radiosonde data will at least yield a not-too-bad representation of the *n* profile in ducting conditions.