Improving the bias characteristics of the ROPP refractivity and bending angle operators

The bending angle observation operator (forward model) currently used to assimilate radio occultation (RO) data at the Met Office, the European Centre for MediumRange Weather Forecasts (ECMWF) and other centres is the same as is included in the Radio Occultation Processing Package (ROPP), along with the corresponding tangentlinear and adjoint code. The functionality of this package will be described in another paper in this issue. The mean bending angle innovations produced with this operator using Met Office background fields show a bias that oscillates with height and whose magnitude peaks between the model levels. These oscillations have been attributed to shortcomings in the assumption of exponentially varying refractivity between model levels. This is used directly in the refractivity operator, and indirectly to produce forward-modelled bending angles via the Abel transform. When the spacing between the model levels is small, this assumption is acceptable, but at stratospheric heights where the model level spacing is large, these biases can be significant, and can potentially degrade analyses. This paper provides physically based improvements to the functional form of refractivity with height. These new assumptions considerably improve the oscillatory bias, and a number of approaches for practical implementation of the bending angle operator are provided.


Introduction
A key feature of radio occultation (RO) data is that the raw observations of excess phase should be unbiased, due to the use of an atomic clock onboard the low Earth orbiting (LEO) RO receiver.These raw measurements, however, are not straightforward to assimilate into a numerical weather prediction (NWP) system, and the raw data are usually preprocessed into bending angles or refractivities, which are then disseminated on the Global Telecommunication System (GTS).
Data assimilation (DA) is the process of producing a statistically optimal "analysis" which is used as an input to an NWP forecast system.The assimilation step blends information from observations and short range (e.g. 6 h) forecast fields, i.e. the background.Mathematically, the basic, timeindependent DA problem is defined as finding the value of x which minimises the following cost function, J : where x is the model state vector, x b is the background state vector (i.e. the "first guess"), y is the observation vector, H is the non-linear observation operator, also called the "forward model", (in a 4D-Var system the operator would include integration of the forecast model) and B and R are the background and observation error covariance matrices, respectively.The first term is evaluated in model space, and the second term is evaluated in observation space.The observation operator is the calculation of the simulated observation which would be measured given the atmospheric state of a model field.Satellites measure quantities such as radiance, excess phase, and not simply atmospheric state quantities such as temperature or humidity.So, the forward model may be fairly complex, even if the observations are pre-processed into quantities more closely related to the model state variables.It should be emphasised that in variational DA, the cost function depends on the "innovations", i.e. y − H (x) and not simply the observations themselves.For this reason, it is just as important to ensure that the forward model H is accurate as it is to ensure that the observations are of good quality.This paper will discuss improvements to the RO refractivity and bending angle forward models which are used at several NWP centres and form part of the Radio Occultation Processing Package (ROPP), henceforth referred to as the "ROPP operator" for brevity.
In the case of refractivity assimilation, the observed refractivity values describe the atmosphere at particular points (with some degree of spatial correlation), so interpolation of refractivity to this point is necessary.Forward-modelled bending angles, however, depend on the entire model atmosphere above the tangent point (Fjeldbo et al., 1971).For this reason, the variation of refractivity with height needs to be known at all heights, from the tangent point to the top of the model, including all points between the model levels.This is necessary in order for the Abel integral, which calculates bending angles from refractivity, to be evaluated (see Sect. 3).
The ROPP operator is based on Healy and Thepaut (2006).This assumes exponentially varying refractivity, N, as a function of x between model levels i and i + 1.The independent variable x is the product of the refractive index n and the distance from the local centre of curvature of the Earth r, i.e. nr: where This ensures continuity at the model levels: With some further approximations, this variation of N with height can allow the bending angle to be calculated via the Abel transform, resulting in a difference of error functions, see Eq. ( 9) (Healy and Thepaut, 2006).To a first approximation the exponential assumption seems reasonable as the refractivity is given by where P is pressure, T is temperature, P w is the partial pressure of water vapour and c 1 and c 2 are empirical constants (Smith and Weintraub, 1953).In a dry atmosphere, the first term in Eq. ( 4) effectively represents the mass field.Where the temperature is constant, the hydrostatic equation dP /dz = −ρg implies that P falls exponentially with height: P (z) = P i exp (−gz/RT ).Therefore, N also falls exponentially.This behaviour can be seen in Fig. 1 for a typical example of how model refractivity varies with height.If the spacing between model levels is sufficiently small, this assumption can produce reasonable refractivities between the levels, and hence bending angles.If the spacing is large, however, the innovation statistics show features which indicate failings of this assumption.This paper aims to address refinements to the form of the refractivity with height used in the refractivity and bending angle forward models.

Refractivity
Currently, bending angles are assimilated operationally at the Met Office but from 2006 to 2010 refractivity data were assimilated, and some NWP centres continue to assimilate refractivity operationally.To forward model refractivity at an observation height which lies between two model levels, the same exponential assumption was applied (i.e.Eq. 2, but in terms of geopotential height).This is equivalent to performing a linear interpolation of ln (N) between the two model levels surrounding the observation height. where With this assumption applied, the innovation statistics ((O − B) /B) are plotted in Fig. 2. All plots in this paper have had Met Office quality control applied to reject potentially poor quality observations (Rennie, 2010).Note that in The refractivity between model levels is calculated using Eq. ( 2).Typical heights of the model levels are overlaid as horizontal lines.
this context, B denotes the simulated observations forwardmodelled from backgrounds, H (x), on observation levels and not the background error covariance as above.
The values of (O − B)/B are calculated for each observed profile (i.e. in observation space), where the background profiles are horizontally interpolated from full-resolution (70level) Met Office fields.The (O −B)/B values are then vertically interpolated linearly onto a fixed grid with 100 m spacing for statistics to be calculated (mean and standard deviation), thus allowing profiles with different sets of impact heights to be included in the statistics.All plotted (O −B)/B statistics in this paper are calculated this way.
The bias above ∼ 45 km should be ignored as it is due to a Met Office-specific model temperature bias, which is anticipated to improve with an upcoming model upgrade.Similarly, the growing negative bias above ∼ 17 km relates, at least partly, to a bias arising from the handling of Met Office levels in the refractivity forward model.This broad bias is potentially problematic, but is specific to the Met Office.The cause is understood and is being addressed but is largely independent of the main topic of this paper, so will be ignored to avoid complicating the discussion.
The general issue that will be addressed here is the smallscale undulation that is present in the bias and is most noticeable between 25 and 45 km.The origin of these fluctuations is clear when the model levels are overlaid, as in Fig. 2.
It can be seen that the magnitude of the oscillatory signal is smallest when the observations are close to the model levels and largest in between.This is a real bias and not a feature of the plotting (the plotting routines have no knowledge of the heights of the model levels, and work entirely in observation Figure 3. Bending angle innovations from the same period as Fig. 2, with typical model levels overlaid.The functional form of refractivity used in the Abel integral is Eq. ( 2).Note that the model levels are plotted on geopotential heights and not converted to impact heights.space).In a DA system (3D-Var for simplicity), the cost function takes the form of Eq. (1).Therefore, the oscillations in the innovations (y − H (x)) will be present in this quantity, and hence they will introduce biases into the DA system.
The origin of these oscillations is apparently the exponential assumption between model levels.

Bending angle
The bending angle forward model is much more sensitive to subtle changes in the model background and the form of dN (x)/dx which is integrated above the tangent height, i.e. the bending angle depends on the vertical gradient of the refractivity.Therefore, it is no surprise that the bending angle statistics show the oscillatory bias even more strongly in Fig. 3.
These statistics are calculated in a similar way to refractivity (above), but the values of (O − B)/B for each profile are interpolated to a fixed grid of impact heights (impact parameter minus the local radius of curvature) rather than geopotential heights.These fixed heights are spaced by 100 m.Plotting bias statistics with coarse vertical binning (e.g. 1 km) can hide these features, so we encourage other NWP centres to follow this methodology to avoid overlooking similar oscillations.
The bending angle as a function of impact parameter α(a) is given by the Abel integral (Fjeldbo et al., 1971;Melbourne et al., 1994;Kursinski et al., 1997): where x = nr as before.By assuming exponential refractivity, and assuming x − a, the bending angle contribution from a single layer is given by Healy and Thepaut (2006): The implementation of the error function uses an accurate fit (Eq. 7.1.25, Abramowitz and Stegun, 1965) to minimise the computational cost.
Because this is an integral from the tangent height upwards and is weighted most strongly close to the tangent point by the denominator, it is expected that if the assumption of exponential refractivity between model levels is less than ideal, then the magnitude of the bias would be smallest close to the model levels, where N(x) is the best representation of the model field (i.e.without any additional distortion from the vertical interpolation), and largest in between.This can be seen in Fig. 3, though unlike the refractivity statistics, the oscillations in the bending angle bias are not symmetric about the centres of the layers.Interestingly, these oscillations do not appear so prominently in the equivalent statistics from the European Centre for Medium-Range Weather Forecasts (ECMWF) -results suggest that this is due to the higher vertical resolution (more than two times larger) in the stratosphere compared to the Met Office, making the exponential assumption more accurate.This is discussed later in this paper.
The bending angle operator proposed by Cucurull et al. (2013) assumes a cubic representation of refractivity as a function of height.This implementation ensures that the vertical refractivity gradients are continuous.The Abel integral is then computed using the trapezoidal rule.In our tests (results not presented here), the oscillatory biases in the innovations were increased for both refractivity and bending angles using this form of N(x), though in our tests the Abel integral was solved analytically rather than numerically.
We therefore seek a new form of refractivity with height as an improvement to the exponential assumption.This can be applied in a number of ways: -Use a more physical function of N(x) as the best approximation, or "reference", between model levels and integrate this or an approximation to it.
-Apply a simple polynomial correction term to the exponential to bring it closer to the reference.
-Use "pseudo-levels"; i.e. evaluate the reference on hypothetical intermediate levels and apply the existing exponential assumption to integrate between these model levels.
These methods all require a best guess for N(x).This should preferably satisfy the following criteria: -N (x) should be continuous at model levels.
-It should have a physical basis.
-It should take information from as few model levels as possible.
-It should include atmospheric moisture.
-It should not be prohibitively costly.

Improved form of N (z)
A form of N (z) used by Healy and Eyre (2000) assumes exponentially varying specific humidity, linearly varying temperature and hydrostatic pressure.This form will be considered as the best guess, or "reference" refractivity between model levels in this paper.In the troposphere, where moisture is most prevalent, the model levels are close together, so the exact form of humidity variation with height is not critical, but exponential variation usually produces a more realistic humidity profile than linear variation in individual cases.The original paper used linear variation of the virtual temperature to obtain the hydrostatic pressure.Here, we use the temperature itself as even at the surface, the difference is rarely more than 1 % and rapidly decreases with height, so in the upper-troposphere-lower-stratosphere, the differences will be negligible.The specific humidity is, however, used in the moist term of the refractivity equation (note that the virtual temperature should be used to compute the geopotential heights on pressure coordinates): where is the ratio of the molecular mass of water vapour and dry air and c 1 and c 2 are as in Eq. ( 4).The reference specific humidity (Q), temperature (T ) and pressure (P ) are defined to behave as .
Between model levels i and i + 1, η i is the inverse scale height of the humidity, β i is the vertical gradient of temperature within the layer, g is the gravitational acceleration and R is the gas constant for dry air.Note that this form of P (z) is different from what is assumed in a previous stage in the Met Office forward model for refractivity; in order to get all model variables on one set of the staggered levels, the Exner pressure values, = (P /P 0 ) R/c p , are interpolated linearly from their native levels.This discrepancy results in the hydrostatic integral producing a discontinuity in N at the model levels.A solution is to replace the temperature gradient β in the expression for pressure (Eq.11), with a value σ that enforces continuity, i.e. where A slightly different version of the continuity correction utilises a factor which scales the pressure linearly within the layer to force continuity.The computed refractivities are almost identical for the two methods, but we choose to proceed with the neater σ correction in this description, as this is the formulation that will form part of the ROPP package.At the Met Office, the alternative formulation is likely to be followed operationally for flexibility, though we emphasise that the underlying assumptions are consistent between these approaches, i.e. the same reference refractivity variation is being approximated.
As stated above, if the forward model handles the model variables consistently throughout, this correction term should not be required.When Eq. ( 10) is used in the refractivity forward model, the vertical profile of the bias becomes significantly smoother, though a small oscillatory signal remains, albeit with opposite curvature at 30 to 40 km.See Fig. 4.
For bending angles, the independent variable is x = nr = n(r curv + z).Because the refractive index is close to unity even near the surface (where n 1.0003), the variation of the refractivity between model levels can reasonably be written in terms of z − z i or x − x i interchangeably.Also, the change to the vertical refractivity gradient arising from this change of variable has been investigated in computations for a small number of cases and the differences are very small.Interchanging these independent variables is only reasonable if nr is monotonic, which is ensured by rejecting observations below any model levels for which the model nr decreases with height.
This approach satisfies the criteria specified in the introduction to Sect. 3.Although we specify that N (z) must be continuous, this new approach does not ensure continuity of dN/dx, which is the quantity integrated in the Abel transform.The importance of this is thought to be small relative to the biases caused by the exponential assumption, and Appendix B contains a specific example and a general demonstration that as long as N is continuous at the model levels, the resulting bending angle profile will also be continuous, regardless of the continuity of dN/dx.

Practical considerations
Two situations can arise where the calculated refractivity is undefined.The first involves the humidity inverse scale height η, defined as In the Met Office 4D-Var system, negative specific humidities can occur throughout the minimisation.This will clearly cause an undefined value of η i , and hence N (z).This is avoided by assuming that Q(z) varies linearly within the layer should the humidity at one of the surrounding levels be negative.In the ROPP package, a positive minimum value of specific humidity is enforced (10 −6 kg kg −1 ).
The second situation is when the temperatures are identical at each of the surrounding levels.In this isothermal case, we initially consider Eq. ( 11).This means that β = 0 and hence P (z) is indeterminate.In this case we therefore replace the expression for P (z) in Eq. ( 11) with its limit as β → 0, namely Knowing that in a dry, isothermal atmosphere the pressure varies exponentially in accordance with the hydrostatic equation, we ensure continuity by replacing the inverse scale height as follows:

Options
Three possible approaches to implement an improved bending angle operator based on the hydrostatic form of the refractivity are presented here.These approaches each have advantages and limitations, and the choice of approach to be implemented will depend on the particular application, including restraints on computational cost.

Expansion of N(x)
If we assume a dry atmosphere, the refractivity reduces to (in terms of x): The hydrostatic pressure will not necessarily be continuous between model levels, so a σ (Eq.13) replaces β in the exponent to preserve continuity of P and hence N: This can be expanded in powers of (x − x i ) to give a correction factor to the exponential: This functional form can also be obtained if instead it is assumed that k i varies linearly within the layer.These two approaches, including the calculation of A and B are described in detail in the Appendix, and their resulting innovation statistics are almost identical.If the moist term is added, this form, i.e.Eq. ( 17), cannot be easily obtained.To use this dry form, a cut-off height is needed (e.g. 12 km), below which, an approach is used that does not require the assumption of a dry atmosphere, such as the existing exponential variation.At these heights, this assumption is reasonable as the model levels are more closely spaced.
The innovation statistics using Eq. ( 19) and the coefficients from the second approach described in Appendix A up to the quadratic term in the series are shown (with no cut-off applied) in Fig. 5.The oscillations in the mean innovations are reduced considerably compared to Fig. 3.There is still an oscillatory feature present in the bias, but now the magnitude is greatest close to the model levels.This may be due to discontinuities in the refractivity gradient, though this has not been investigated.

Polynomial correction
The exponential form of N(z) can be modified by additional terms to better approximate the "reference" refractivity, including the moist term.For example (redefining A i and B i ): Bending angle innovation statistics using an integrable approximation to the dry hydrostatic refractivity at all heights, i.e.Eq. ( 19).See Sect.A3 for full details.
This could be used to give a very good approximation to the "reference" (if we know it), and can easily be integrated in the Abel transform, resulting in extra terms in addition to the error function.Figure 6 shows typical differences between the hydrostatic refractivity (Eq.10) and the exponentially varying refractivity between two model levels, as well as a quadratic approximation to this difference as described below.As a polynomial correction is a fit to the difference between the "reference" (i.e. the hydrostatic refractivity) and the exponential form, this difference must be specified at a number of points that is commensurate with the degree of the correction in order to fully determine the fit.For the quadratic example shown in Fig. 7, the values of the quadratic correction at the two surrounding model levels are set to zero to ensure continuity, and the difference between the corrected hydrostatic and exponential forms at the centre of the layer (i.e. the horizontal dotted line) is used to provide the remaining information to fully determine the quadratic correction.
For continuity at z i+1 , the following relation must hold, since k i is still given by Eq. ( 3): The value of the quadratic at its turning point is set to be equal to the difference between the hydrostatic and exponential forms of the refractivity at the layer midpoint.This is reasonable to assume, as from visual inspection the differences are approximately quadratic (Fig. 6), and hence fairly symmetric about the midpoint.The turning point of the correction is found by setting the first derivative of the correction to zero: If the turning point is close to the middle of the layer, we can substitute Eq. ( 22) into the expression for the quadratic Figure 6.The difference of the corrected hydrostatic refractivity and the exponentially varying refractivity between a single pair of Met Office model levels (horizontal lines).Also shown is a quadratic approximation to this difference, as described in the text.The horizontal dotted line shows the midpoint of the layer where the quadratic correction is set to be equal to the difference of the hydrostatic and exponential values.correction at the midpoint, where N hyd_mid and N exp_mid are the refractivity values at the middle of the layer calculated using the hydrostatic and exponential approaches, respectively.Substituting A i from Eq. ( 21), we obtain a value for B i : Inserting this form into the Abel integral results in an additional term in the expression for bending angle (having swapped z − z i for x − x i in an intermediate step): This has been extended to include a cubic term to account for the small asymmetry in N hyd − N exp at the midlayer point.This does not show a significant improvement and leads to a more complicated form of the integral, so the results are not presented here.
The polynomial correction has the advantage that the humidity is accounted for, and the first order behaviour is already accounted for by the exponential, so other reference refractivities could be used to provide updates to the coefficients in the future.Bending angle innovation statistics, using a quadratic adjustment to the exponential form of refractivity with height (Eq.20) to produce a better approximation to the hydrostatic form.

Pseudo-levels
If the "reference" refractivity, including the moist term, is evaluated at intermediate "pseudo-levels" which lie between the model levels (having first calculated Eq. ( 11) on these pseudo-levels, ensuring continuity of the pressure), then the exponential assumption can be accurately applied between these levels (if there are sufficient additional levels), so the current (exponential) operator can simply be invoked multiple times within each model layer.For future changes, this is a flexible approach as the computation only needs to evaluate the refractivities on the pseudo-levels and the Abel integral remains unchanged, hence additional assumptions/simplifications can be avoided and a more sophisticated form could potentially be used.The number of pseudo-levels must be chosen to provide a balance between accuracy and computational cost.It has been found that using just one additional pseudo-level in the middle of each layer gives a good improvement for the associated cost.Two or more equally spaced pseudo-levels only provide very small improvements to the innovation statistics for the single pseudo-level case, so results with just one pseudo-level are presented here.For the layer in which the tangent point lies, the refractivity expression, Eq. ( 10), is used to evaluate N at the tangent height, and at an additional pseudo-level halfway between the tangent point and the next highest model level.The resulting innovation statistics are shown in Fig. 8.
A further use of this method has been to examine the effect of "doubling" the number of model levels by introducing mid-layer pseudo-levels.This is similar to what is described above, but the treatment of the layer in which the tangent point lies is slightly different -the pseudo-level in this layer is at the layer's midpoint, rather than halfway between the Figure 8. Bending angle innovation statistics, using hydrostatic refractivity (Eq. 10, including the moisture term) evaluated on one additional pseudo-level per model layer, and using an exponential function of refractivity with height to evaluate the Abel integral between the model/pseudo-levels.
tangent height and the next model level as was described above.The motivation for investigating this is to explain why the innovations from the L91 ECMWF system (ECMWF, 2007) do not show these oscillations as strongly as in the L70 Met Office (Davies et al., 2005) statistics.At a height of 35 km, where the oscillations in the bias are prominent, the level spacing of the L91 ECMWF model is ∼ 1.5 km, and at the Met Office (L70) it is ∼ 2.9 km, i.e. a factor of ∼ 2 different.Figure 9 shows the innovations when pseudo-levels are used in this configuration.
By comparing Figs. 9 and 3, it can be seen that by doubling the effective number of levels, the oscillations are reduced, and hence this provides an explanation as to why the ECMWF statistics do not display these features as strongly.In other words, the exponential assumption is more acceptable with the L91 resolution, but less so for L70.
Similarly, when the ECMWF levels are thinned by a factor of two, the innovation statistics show the oscillatory bias much more strongly, and is very similar to the Met Office bias structure.This is shown in Figs. 10 and 11.The ECMWF implementation used in these plots is described in Appendix A3 and uses a 12 km cut-off, below which the original operator is used.
Another contributing factor to the smaller oscillatory bias using ECMWF profiles is that the ECMWF height levels are more variable in this region than the Met Office levels and this could lead to the smoothing out of the oscillatory signal, but this effect has not been investigated here.
For reasons of longer-term flexibility and maintenance, this approach is due to be implemented at the Met Office in 2014, whereas the expansion of the dry refractivity (described in detail in Sect.A3) will be implemented in ROPP, though both approaches are based on the same underlying principles.

Conclusions
It has been demonstrated that when the vertical model level spacing is large, the assumption of exponentially varying refractivity leads to systematic negative biases in forwardmodelled stratospheric refractivities and bending angles for which the magnitudes are largest when the observation height lies between the model levels.The use of a more physical form of refractivity as a function of height has been investigated.This function assumes exponentially varying humidity, linearly varying temperature and hydrostatic pressure.Using this function, the magnitude of the oscillatory bias has been reduced considerably in both refractivity and bending angle statistics using Met Office background profiles.Three approaches to implement such an improvement have been suggested: 1. integrate an approximation to the dry-hydrostatic refractivity analytically above a point where the moist refractivity term is negligible; 2. apply a polynomial correction to the exponential to make it a better approximation to the hydrostatic form; 3. evaluate the hydrostatic refractivity on mid-layer pseudo-levels and use the exponential function in the Abel integral between the model/pseudo-levels.
These methods each have their own merits, and these have been stated in the text.
In Appendix A, two methods of approximating the dry hydrostatic form are given and the resulting bending angle statistics are consistent.
The results presented here should provide an improvement to operational DA systems.Usually, RO data is assimilated without a bias correction, and hence acts as an anchor (Poli et al., 2010;Healy, 2008) to correct biased radiance observations.It is anticipated that the reduction of this forwardmodel bias will improve analyses both directly and indirectly via bias correction schemes.Findings reported here could also be used in 1D-Var retrieval chains to improve the quality of the retrieved quantities, as well as reanalysis and climate model validation.

C. P. Burrows et al.: Radio occultation forward models
so that a sudden change in dT / dx would cause a jump in dN/ dx.Note that we assume that N itself is continuous everywhere, and that dN/ dx is finite everywhere.

B1 In particular
To be specific, we assume where N 0 = N (x 0 ) and k 1 > k 0 (> 0) for definiteness (corresponding to a more positive dT / dx above x 0 in the "tropopause" model above).This implies The key point is that the bending angle is continuous at x 0 : α(x + 0 ) = α(x − 0 ) = √ 2πx 0 k 1 N 0 .A secondary point is that the same cannot be said for the derivative of α -indeed, dα/ da (x − 0 ) is formally infinite.In fact, for a just below x 0 , Eq. (B6) implies α(a)−α(x 0 ) = 2 2x 0 (x 0 − a)N 0 (k 0 −k 1 )+O(x 0 −a).(B7) Note that α(a) < α(x 0 ) when k 1 > k 0 .This is because the (x − a) −1/2 factor in Eq. (B1) means that α(a) is dominated by the contribution from N just below x 0 , which in this case is smaller (in magnitude) than N just above it.
Figure B1 shows N , dN/dx and α for a 15 km "tropopause".k 0 = 0.1 km −1 ; k 1 = 0.2 km −1 .The continuity of α at x 0 = 15 km is clear, as is its cusp just below.The refractivity at the "tropopause" is 45 N-units, and the radius of curvature used in the bending angle calculation is 6350 km.

B2 In general
More generally, suppose that there is a jump in dN/ dx at x 0 .Is the bending angle continuous there?

Figure 1 .
Figure1.Logarithm (base 10) of a randomly selected (but fairly typical) vertical refractivity profile, forward modelled from a 70level Met Office background profile.This highlights the approximately exponential behaviour of refractivity with height.The dry and wet terms have also been plotted to show their relative contributions.

Figure 2 .
Figure2.Refractivity innovations from 25 Met Office (6-hourly) model cycles, with observation data from all available RO instruments.The period started with the 00 Z analysis on 1 January 2014.The refractivity between model levels is calculated using Eq.(2).Typical heights of the model levels are overlaid as horizontal lines.

Figure 4 .
Figure 4. Refractivity innovations using the hydrostatic refractivity expression between model levels (Eqs.10 and 11) with typical heights of the model levels overlaid.

Figure 7 .
Figure7.Bending angle innovation statistics, using a quadratic adjustment to the exponential form of refractivity with height (Eq.20) to produce a better approximation to the hydrostatic form.

Figure 9 .Figure 10 .
Figure9.Bending angle innovation statistics, using hydrostatic refractivity (Eq. 10, including the moisture term) evaluated on a doubled-resolution vertical grid and using an exponential function of refractivity with height to evaluate the Abel integral between the model/pseudo-levels

Figure 11 .
Figure 11.Bending angle innovation statistics as per Fig. 10, but with the background profiles thinned to half the vertical resolution of the 91-level ECMWF model.

Figure B1 .
Figure B1.Example profiles of (from left to right) N, −dN/ dx and α when there is a discontinuity in dN/ dx at x 0 = 15 km.(For this plot, N = 10 6 (n − 1) as usual.