8.  Atmospheric Delay Models

8.1  Description of the atmospheric  delay

As the GPS signal travels from the satellite to the receiver, it propagates through the atmosphere of the earth, where it is retarded and its path changed from a straight line to a curved one.  If we take the simplified mathematical model for the observable to be one in which the signal is assumed to be propagating in a straight line and at the speed of light in vacuum, then the "atmospheric delay" is defined to be the difference between the true electrical path length and this assumed straight-line length. Using this definition, the atmospheric delay is a term to be added to the simplified model. 

The atmospheric delay in the zenith (i.e., vertical) direction varies from about 6 to 8 nanoseconds (190 to 240 cm, or 10-12 cycles of phase at L1-band) depending on meteorological conditions and site location.  The atmospheric delay increases with decreasing elevation angle approximately with the cosecant of the elevation angle, so that the atmospheric delay at an elevation angle of 20 degrees may be from 30-36 cycles of L1 phase.  

The atmospheric delay is usually broken down into two components.  The first component is due to the mixture of all constituents, but it is assumed that the mean molar mass of these constituents is equal to the mean molar mass of only the " dry" (all except water vapor) constituents.  Assuming that the atmosphere is in hydrostatic equilibrium, the "zenith delay" due to these components is very well modeled (standard deviation of approximately 0.5 mm) using the surface pressure, which represents the total weight of the atmosphere.  This component of the atmospheric propagation delay is usually termed the "dry" or "hydrostatic" delay, and accounts for nearly all ( 90-100%) of the atmospheric propagation delay.  

The second component of the atmospheric delay is due to water vapor, and includes a correction for the "dry mean molar mass" assumption used to derive the dry delay (see above).  This component of the atmospheric propagation delay is called the "wet delay" and is equal to zero if there is no water vapor present anywhere along the path of the signal.  However, there usually is  water vapor present along the path of the signal and it is poorly predicted using measurements of conditions at the site alone.  This difficulty is caused by the "unmixed" condition of atmospheric water vapor, which means that the water vapor is present in "blobs" throughout the troposphere.  Because of this condition, models for the wet delay are notoriously inaccurate and can have RMS errors of several cm (zenith), out of a total (zenith) wet delay of 0-40 cm. 

8.2  Algorithms for the atmospheric propagation  delay

The atmospheric propagation delay is the implemented in the following manner: 

    	ATDEL(EL) = DRYZEN * DRYMAP(EL) + WETZEN * WETMAP(EL) 
  
where EL is the elevation angle of the satellite, DRYZEN is the dry zenith delay, WETZEN is the wet zenith delay, DRYMAP is the "mapping function" for the dry delay (see below) and WETMAP is the mapping function for the wet delay.  A mapping function is a mathematical model for the elevation dependence of the respective delays.  The mapping functions (for both the dry and the wet terms) are approximately equal to the cosecant of elevation, but there are significant deviations from this "cosecant law" due both to the curvature of the earth and the curvature of the path of the GPS signal propagating through the atmosphere. 

Many expressions for the four terms DRYZEN, DRYMAP, WETZEN, and WETMAP appear in the scientific literature.  For microwave observations there is little controversy about the expressions for the dry zenith delay.  Since the wet zenith delay cannot be accurately modeled from surface measurements, the expressions used are not critical.  The choice of mapping function can be important, however, particularly for elevation cutoff angles below 15 degrees.  The model to be used in GAMIT for each of the four terms is specified by keywords in the sittbl. (see Section 5.2), a portion of which is shown below with the default choices and values of the surface pressure, temperature, and relative humidity.

SITE              FIX  --COORD.CONSTR.-  DZEN  DMAP  WMAP  ---MET. VALUE----  WFILE
CATO Castro Peak  NNN  100.  100.  100.  SAAS  SAAS  CFA  1013.25  20.0 50.0  NONE

DZEN:
 
The default model [ SAAS ] for the dry zenith delay is that described by Saastamoinen [1972].  If there is a W-file of meterological data, the actual pressure at the station, as interpolated from the file is used to determine the delay.  If there is no W-file ( NONE, the default) the model takes the surface pressure given in the sittbl. and reduces it to pressure at the station using the station's height and small corrections for gravitational variations with latitude.  You can also input directly a constant value for the dry zenith delay by specifying the model as CNNN where NNN is the range-equivalent of the delay in centimeters; e.g., C230 sets the delay as 230 cm.  This feature is useful mainly for error analysis.

WZEN:
   
The default model [ SAAS ] for the wet delay is also from Saastamoinen [1972].  It makes use of both the surface temperature and the relative humidity/dewpoint, either from the meteorological data file or from the default values.  The default values of temperature and humidity (20 degrees C and 50%) serve to make the a priori delay "reasonable" for temperate climates.  As for the dry zenith delay, you may set a constant value using CNNN as the input model except that the constant delay is given in millimeters (NOT centimeters).  This feature is useful because, as mentioned above, the true wet zenith delay is often not correlated with surface measurements.  You can use C000 (i.e., zero wet zenith delay) and estimate a zenith delay correction which, if the correct pressure was used for the dry zenith delay, will represent the mean wet zenith delay.  An additional option for the wet delay is to determine the value from water-vapor radiometer data, as described in Section 8.5 below; in this case set the WZEN code to WVR.  

DMAP:

The default mapping function for the dry delay is CfA-2.2 [ CFA] , developed at the Harvard-Smithsonian Center for Astrophysics [Davis et al., 1985] for analysis of VLBI data.  The original model uses input values of tropospheric temperature lapse rate and tropopause altitude as well as surface pressure, temperature, and vapor pressure, but GAMIT assumes nominal values for the first two quantities.  New mapping functions have recently been developed that represent measurable improvement over CfA-2.2 for elevation angles below 10 degrees.  The version developed by Arthur Niell of MIT's Haystack Observatory  [Niell, 1996] uses standard values of meteorological parameters and has been recently added to the GAMIT code.  The hydrostatic (dry) formulas may be invoked by specifying NMFH for DMAP.  (The other recent mapping function is the MTT model developed by Herring [1992], but this has not yet been coded in GAMIT.)  GAMIT also supports two older dry mapping functions which were commonly used for analysis of microwave observations prior to 1985—one developed by Marini and Murray [1974] [MARI ] and one by Chao [1972] [ CHAO ]

WMAP:  

'CFA ' is the default, but the Niell model [ NMFW ] is superior and Marini [ MARI ] or Chao   [ CHAO ] models may be used.  

Formulas coded

  P - Total surface pressure, mbars 
  T - Surface temperature, Kelvins 
   e - Surface water vapor pressure, mbars 

<<<SAAS>>> 

	çf = site latitude 
	H   = ellipsoidal height, km  
	F(f,H) = 1 - 0.00266 * cos(2 f ) - 0.00028 * H 
	SAAS(dry) = 0.2277 * P / F(f,H) in centimeters 
	SAAS(wet) = 0.2277 * e * (0.05 + 1255 / T) / F(f,H)   

 <<<MARI>>>
  
	BETA  = 2.644E-03 * exp(-0.14372 * H) / F(f,H)
	KAPPA = BETA / [SAAS(dry) + SAAS(wet)] 
	EL    = Elevation of satellite 
	MARI = (1 + KAPPA) / { sin(EL) + [KAPPA/(1+KAPPA)]/[1+SIN(EL)] } 

<<CFA >>> 

	A = 0.001185 * [1 + 0.6071E-04 * (P - 1000)
	                               - 0.1471E-03 *  e 
	                              + 0.3072E-02 * (T - 293)  ]  

	B = 0.001144 * [1 + 0.1164E-04 * (P - 1000)
                                       + 0.2795E-03 *  e
                                       + 0.3109E-02 * (T - 293)  ] 
 
	C = -0.0090   CFA = 1/ { sin(EL) + A / [ tan(EL) + B / ( sin(EL) + C ) ] } 

8.3  Estimating a zenith delay parameter

Since the water vapor contribution to atmospheric delay is poorly modeled using surface meteorological data, GAMIT allows estimation of corrections to the zenith delay.  The partial derivative of phase or pseudorange with respect to the zenith delay parameter is simply the mapping function, approximately equal to the cosecant of the elevation angle of the satellite as viewed from the station.   For stations in a regional network the elevation angles viewing a particular satellite will be nearly equal, producing high correlations among the estimated zenith delays.  Even for stations a few meters apart, however, separate zenith delays can be estimated without causing numerical problems in SOLVE; although the uncertainties in all of the zenith-delay parameters will be large, the relative values of the estimates themselves can be trusted.  To extract from the solution the uncertainties of the differences, however, you should fix or tightly constrain one of the values. 

The model for zenith delay can take the form of a single parameter for each station and session, or a piecewise linear function of zenith-delay over the session.  In the latter case, the tabular points of the function can be constrained using a first-order Gauss-Markov process.  Controls for estimation of zenith delay are input via the sestbl. and/or sittbl.  The sestbl. inputs, adopted as common to all stations are as follows:

Zenith Delay Estimation = YES  ; YES/NO
Number Zen = 4                 ; number of zenith-delay parameters
Zenith Model = PWL             ; PWL (piecewise linear)/CON (step)
Zenith Constraints = 0.50      ; zenith-delay a priori constraint in meters
Zenith Variation = 0.02 100.   ; zenith-delay varation, tau; units m/sqrt(hr), hrs

Specifying Number Zen = 1, and either PWL or CON for Zenith Model will invoke a single parameter for the zenith delay over the session.  The best representation of zenith-delay variations is usually accomplished with a new zenith delay parameter for every 2–6 hours during the day (the gain from 6 hrs to 2 hrs is usually small and will increase running time for SOLVE considerably); thus with an 8-hr observation span, you might set Number Zen=3 to get two 4-hr segments (the PWL representation includes tabular points on each end of the observation span), or with a 24-hr span, set Number Zen=9 to get 3-hr segments. The overall zenith constraint should be set loose enough to encompass comfortably any error in the wet delay; 0.5 meters is the default and reasonable.  The variation is specified as parameters of a first-order Gauss-Markov process.  The first value in Zenith Variation is the point-to-point variation allowed, in units of meters.  The second value is the correlation time ( tau ) in hours.  Setting tau long compared to the observation span results in a random walk process, which is both reasonable and easy to interpret (and has the practical advantage of persistence with large error bars for spans with few observations). the easiest to understand   The default value of 100. hrs accomplishes this for 24-hr spans.  Setting tau equal to a value short compared with the tabular point interval will result in a white-noise process for the variation in tabular points; in this case the contraint will be applied with respect to the default model value rather than the value of the last tabular point. 

There is an additional entry in the sestbl.,

Tropospheric Constraints = NO   ; YES/NO  

which invokes a spatial constraint on the zenith-delay parameters.  This constraint can be useful for tying together the zenith-delay adjustments for closely-spaced sites in a network. This feature was coded originally for a single zenith delay, however, and does not yet work for time-dependent models. 

Different values for the zenith constraints can be invoked using sittbl. entries, shown below:

SITE              DZEN  DMAP  WMAP  ---MET. VALUE----  WFILE  ZCNSTR  ZENVAR  ZENTAU
CATO Castro Peak  SAAS  SAAS  CFA  1013.25  20.0 50.0  NONE   0.5000  0.005   100. 
LOVE Loma Verde   SAAS  SAAS  CFA  1013.25  20.0 50.0  NONE   0.5000  0.005   100. 
PVER PVER7268     SAAS  SAAS  CFA  1013.25  20.0 50.0  NONE   0.5000  0.020   100. 
SAFE San Fernando SAAS  SAAS  CFA  1013.25  20.0 50.0  NONE   0.5000  0.020   100. 
YKNF YELLOWKNIFE  SAAS  SAAS  CFA  1013.25  20.0 50.0  NONE   0.5000  0.020   100. 
AROG ARO GPS      SAAS  SAAS  CFA  1013.25  20.0 50.0  NONE   0.5000  0.020   100. 
RICH RICHFRPA     SAAS  SAAS  CFA  1013.25  20.0 50.0  NONE   0.5000  0.020   100. 
TROM TROMGPSM     SAAS  SAAS  CFA  1013.25  20.0 50.0  NONE   0.5000  0.020   100. 

It is currently a requirement that the number of zenith delays in the session be the same for all stations. 

8.4  Estimating gradients

The effects of azimuthal asymmetry in the atmospheric delay are not included in MODEL but may be estimated in SOLVE.  The coded partials imply a model of the form  

ATDEL(EL,AZ) = GRADNS*AZMAP(EL)*COS(AZ) + GRADEW*AZMAP(EL)*SIN(AZ)
	
where EL is the elevation angle, AZ the azimuth, and AZMAP the mapping function for gradients, given by 

	AZMAP = 1./ (SIN(EL)*TAN(EL) + C

and C is a constant equal to 0.003  [Chen and Herring, 1996].  Since the gradient parameters, GRADNS and GRADEW, have small and non-intuitive values near the zenith (i.e., for AZMAP = 1 ), we rescale them to represent the difference between the the north (or east) and south (or west) delay at 10 degrees elevation.  At 10 degrees the rms scatter of gradients observed from VLBI observations are about 5 mm.  Our default a priori constraint is 30 mm.

8.5  Water vapor  radiometer (WVR)  data 

 A water-vapor radiometer is a multichannel radiometer which measures the brightness temperatures of the atmosphere.  The main contributions to this radiation below 100 GHz are atmospheric oxygen, water vapor, and liquid water.  The several frequencies of the WVR are used to solve for these various components and the wet delay.  The algorithm for the delay currently implemented is  

       	Wet delay = A + B * T1 + C * T2  

where A, B, and C are user-defined constants, and T1 and T2 are the brightness temperatures (not linearized) at two frequencies.  The name of the WVR  data file (Z-file) must also be entered.  The Z-file is of the format:  

This is a comment line      
   .     
   .       
This is the last comment line 
END 
Header line #1 
Header line #2 
YYYY DDD HH MM  YYYY DDD HH MM  SITE
Header line #3
Header line #4 
YYYY DDD HH MM SS   AZ EL T1 T2 TS    
    .      
   .     
   .  

The actual entries are free format.  They include the year, day-of-year, hours , minutes and seconds for the WVR observation, the azimuth and elevation (in degrees) of the observation, the brightness temperatures (in Kelvins), and the surface temperatures (in degrees Celsius).  In order to interpolate the brightness temperatures, they must be converted to  quivalent zenith brightness temperatures.  This is done in the following manner:  Saturation effects are taken into account by determining the optical depth t using   

  		Ti = T0 exp (-ti ) + Teff [1 - exp(-ti)]  

where ti is the optical depth at the frequency of the brightness temperature Ti,  T0=2.9 K is the big-bang background temperature, and Teff is the "effective" temperature of the atmosphere [Chandraesker, 1960; Wu, 1979].  We use the approximations  

   		Teff = 0.94 * Ts 
 
where Ts is the surface temperature in Kelvins [Wu, 1979].  The optical depth is then mapped to the zenith using 

    t = t (zen) * WETCHAO(elevation)
  
where the function WETCHAO is the wet Chao mapping function given above.  The zenith brightness temperatures are then calculated from the zenith opacities using the expression above.  These zenith brightness temperatures are then interpolated to the epoch  of the GPS measurement, and the zenith delay is calculated using the user-defined WVR constants.  The wet delay away from the zenith is computed using the Chao wet mapping function.  Note: For consistency with the interpolation scheme, the user must choose 'CHAO' for the wet mapping function if WVR is chosen for the zenith wet delay.  If this is not done, the program will display a warning and the wet mapping function will be changed to 'CHAO'. 

8.6 References  

Chandrasekhar, S., Radiative Transfer, Dover, New York, 1960. 

Chao, C. C., A model for tropospheric calibration from daily surface and radiosonde balloon measurements, Tech. Mem. 391-350, Jet Propulsion Laboratory, Pasadena, California, 1972. 

Davis, J.L., T.A. Herring, I.I. Shapiro, A.E.E. Rogers, G. Elgered, Geodesy by radio interferometry:  Effects of atmospheric modeling errors on estimates of baseline length, Radio Science, 20, 1593-1607, 1985.

Chen, G., and T. A. Herring, Effects of atmospheric azimuth asymmetry in the analysis of space geodetic data, in preparation, 1996.

Herring, T. A., Modeling atmospheric delays in the anlaysis of space geodetic data, Proceedings of the Symposium on Refraction of Transatmospheric Signals in Geodesy, Publications on Geodesy, No. 36, ed. by De Munck and Spoelstra, 247, Netherland Geodetic Commission, Delft, Netherlands, 1992. 

Marini, J.W., C.W. Murray, Correction of radio range tracking data for atmospheric refraction at elevations above 10 degrees, unpublished memorandum, 1974. 

Niell, A. E., Global mapping functions for the atmospheric delay at radio wavelengths, Journal of Geophysical Research, 101, 3227–3246, 1996. 

Treuhaft, R. N., and G. E. Lanyi, The effect of the dynamic wet troposphere on VLBI measurements, Geodesy and Geophysics Preprint No. 141, Jet Propulsion Laboratory, Pasadena, 1986. 
 
Saastamoinen, J., Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites, in The Use of Artificial Satellites for Geodesy, Geophys. Monogr. Ser., Vol. 15, edited by S.W. Henriksen et al., pp. 247-251, American Geophysical Union, Washington, D.C., 1972. 

Wu, S.C., 1979: Optimum frequencies of a passive microwave radiometer for tropos pheric path-length correction.  IEEE Trans. Ant. Prop., AP 20, 181-187.